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Abstract 


This work is about the inverse dynamics of underactuated flexible mechanical systems 
governed by quasi-linear hyperbolic partial differential equations subjected to time-vary- 
ing Dirichlet boundary conditions that are enforced by unknown, spatially disjunct, hence 
non-collocated Neumann boundary conditions. The governing equations are first derived 
abstractly in more detail before various mechanical systems are presented aligning with 
the postulated formulation. Therefore, geometrically exact theories, governing the large 
motion of slender structures such as strings and beams but also general three-dimensional 
continua, will be derived. 


Commonly, initial boundary value problems that occur in non-linear structural dynam- 
ics are solved by applying sequential space-time discretization methods. This approach 
is therefore typically based on a discretization in space by means of finite elements, fol- 
lowed by an appropriate discretization in time mostly based on finite differences. A brief 
survey of this type of sequential space-time integration methods for the initial bound- 
ary value problem at hand, is given first in context of the direct dynamics problem. Le. 
the pure Neumann boundary problem will be considered first, before different possibili- 
ties of imposing Dirichlet boundary conditions in general, will be discussed afterwards. 
Based on this, the inverse dynamics problem will be introduced in context of spatially 
discrete mechanical systems subjected to rheonomic holonomic servo-constraints. A de- 
tailed analysis of this type of constrained systems aims to elaborate the fundamental 
distinctions of servo-constraints to classical contact-constraints. Consequences thereof, 
regarding the construction of numerically stable integration methods, will be addressed 
likewise before numerous selected examples will be given. 


Due to the highly restrictive applicability of solving the inverse dynamics problem se- 
quentially in space and time, an in-depth analysis of the underlying initial boundary 
value problem is being pursued. Especially by exposing the underlying hyperbolic struc- 
ture of the governing partial differential equations, further insights into the problem at 
hand are anticipated. Enlighting the resulting mechanisms of wave propagation within 
continuous structures will pave the way to a numerically stable integration of the in- 
verse dynamics problem. Thus, e.g. a method that is based on the integration of partial 
differential equations along characteristic manifolds is presented. This motivates the de- 
velopment of novel Galerkin methods that can be presented in this work as well. 


These newly established methods will be applied to the feed-forward control problem of 
various mechanical systems. In addition to that, the novel simultaneous space-time inte- 
gration methods are adopted to flexible multibody systems such as e.g. the cooperative 


Abstract 


control of a rigid body through several flexible strings or the control of the end-effector 
of a rotational arm consisting of rigid and flexible arms. 


Selected numerical examples are given underlining the relevance of the proposed simul- 
taneous space-time integration. 
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Kurzfassung 


Diese Arbeit befasst sich mit der inversen Dynamik unteraktuierter, flexibler, mechanis- 
cher Systeme, welche durch quasi-lineare hyperbolische partielle Differentialgleichungen 
beschrieben werden können. Diese Gleichungnen, sind zeitlich veränderlichen Dirich- 
let-Randbedingungen unterworfen, welche durch unbekannte, räumlich disjunkte, also 
nicht kollokierte Neumann-Randbedingungen erzwungen werden. Die zugrundeliegen- 
den Gleichungen werden zunächst abstrakt hergeleitet, bevor verschiedene mechanis- 
che Systeme vorgestellt werden können, die mit der eingangs postulierten Formulierung 
übereinstimmen. Hierzu werden geometrisch exakte Theorien hergeleitet, welche in der 
Lage sind große Bewegungen schlanker Strukturen wie Seile und Balken, aber auch ganz 
allgemein, dreidimensionaler Festkörper zu beschreiben. 


In der Regel werden Anfangs-Randwertprobleme, die in der nichtlinearen Strukturdy- 
namik auftreten, durch Anwendung einer sequentiellen Diskretisierung in Raum und Zeit 
gelöst. Diese Verfahren basieren für gewöhnlich auf einer räumlichen Diskretisierung 
mit finiten Elementen, gefolgt von einer geeigneten zeitlichen Diskretisierung, welche 
meist auf finiten Differenzen beruht. Ein kurzer Überblick über derartige sequentielle 
Integrationsverfahren für das vorliegende Anfangs-Randwertproblem wird zunächst an- 
hand der direkten Formulierung des Problems gegeben werden. D.h. es wird zunächst 
das reine Neumann-Randproblem betrachtet, bevor anschließend ganz allgemein, ver- 
schiedene Möglichkeiten zur Einbindung etwaiger Dirichlet-Randbedingungen diskutiert 
werden. Darauf aufbauend wird das Problem der inversen Dynamik im Kontext räumlich 
diskreter mechanischer Systeme, welche rheonom-holonomen Servo-Bindungen unter- 
liegen, eingeführt. Eine ausführliche Untersuchung dieser Art von gebundenen Syste- 
men soll die grundlegenden Unterschiede zwischen Servo-Bindungen und klassischen 
Kontakt-Bindungen herausarbeiten. Die daraus resultierenden Folgen für die Entwick- 
lung geeigneter numerisch stabiler Integrationsverfahren können dabei ebenfalls ange- 
sprochen werden, bevor zahlreich ausgewählte Beispiele vorgestellt werden können. 


Aufgrund der sehr eingeschränkten Anwendbarkeit der sequentiellen Lösung der inversen 
Dynamik in Raum und Zeit, wird eine eingehende Analyse des vorliegenden Anfangs- 
Randwertproblems unternommen. Vor allem durch die Freilegung der hyperbolischen 
Struktur der zugrundeliegenden partiellen Differentialgleichungen werden sich weitere 
Einblicke in das vorliegende Problem erhofft. Die Erforschung der daraus resultierenden 
Mechanismen der Wellenausbreitung in kontinuierlichen Strukturen öffnet die Tür zur 
Entwicklung numerisch stabiler Integrationsverfahren für die inverse Dynamik. So kann 
unter anderem eine Methode vorgestellt werden, die auf der Integration der partiellen Dif- 
ferentialgleichungen entlang charakteristischer Mannigfaltigkeiten beruht. Dies regt zu 
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Kurzfassung 


der Entwicklung neuartiger Galerkinverfahren an, die ebenfalls in dieser Arbeit vorgestellt 
werden können. 


Diese neu entwickelten Methoden können anschließend auf die Steuerung verschiedener 
mechanischer Systeme angewendet werden. Darüber hinaus können die neuartigen Inte- 
grationsverfahren auch auf flexible Mehrkörpersysteme übertragen werden. Angeführt 
seien hier beispielsweise die kooperative Steuerung eines an mehreren flexiblen Seilen 
aufgehängten starren Körpers oder die Steuerung des Endeffektors eines flexiblen mehr- 
gliedrigen Schwenkarms. 


Ausgewählte numerische Beispiele verdeutlichen die Relevanz der hier vorgeschlagenen, 
in Raum und Zeit simultanen Integration des vorliegenden Anfangs-Randwertproblems. 
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1. Introduction 


How does a flexible body need to be excited on a prescribed boundary, such that the 
motion on a spatially disjunct boundary can partly be specified? This seemingly simple 
question will be pursued in this work and problems, that arise by attempting to solve 
such inverse dynamics problems of spatially continuous systems using classical numeri- 
cal solution strategies, will be pointed out. 


1.1. Motivation and formulation of the problem 

In the sequel, we will focus on infinite-dimensional systems governed by second-order 

quasi-linear hyperbolic partial differential equations of the form 
A-&x-divs(B:9x)=C. (1.1) 


Since we are interested in quasi-linear equations, we explicitly allow the coefficients A, 
B and C to depend on the variables s € Sc Rf andt € T c R,, as well as on the solution 
x: Q m R@ and their first partial derivatives itself, such that 


A,B: Qh R?! and C:ÖH>R where Q=QU {(x, 3x, 3x) : Q m RA. 


For convenience, the space-time domain Q = S x T has been introduced, where S c R” 
represents the spatial and 7 C R, the temporal domain. Note that the spatial dimension 
d > 1, does not necessarily coincide with the dimension of the spatial variable a > 1. 


In context of the classical direct formulation of the problem, the task is to find a solution 
x: Qh R@ satisfying Eq. (1.1), that additionally is subjected to some given Neumann 
boundary conditions on 094 = 05, X T 


Bdsx\|(s,1) car = f(s,t) (1.2) 
and (time-varying) Dirichlet boundary conditions on dQ, = 0S, x T! 


x (scan, = V(S.t) (1.3) 


1 In the special case of spatially one-dimensional systems, i.e. d = 1, we define aS f = {0} and aS, = {L}, for 
some L € R. 
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as well as initial conditions on Qo = S x {0} 
*|(styeoo, =X0 and IX] (s, tJe% = Vo- (1.4) 
Note that for the direct problem the conditions 


Bdsx|(5,1)<02, = 1s, t) (1.5) 


on Qy need to be undetermined, i.e. Dirichlet and Neumann boundaries need to be 
disjunct. 


In contrast to that, the inverse problem seeks to find a function f : Qf > R@ on ðQ rby 
additionally imposing Neumann boundary conditions on dQ, = 0S, x T as given by 


BdsX\(5.1)<a0, = n( 07x, 0;x,x,t) . (1.6) 


Note that the function 7 in general might be given in differential form. 


Whereas in the direct problem given Dirichlet boundary conditions (1.3) on dQ, are en- 
forced by the unknown function ņ : dQ, > Rf, acting on the same boundary 99, the 
inverse problem seeks to find the unknown function f : Qf —> R@ on ðQ f enforcing 
(1.3). This unknown Neumann-boundary conditions can then be considered as Lagrange 
multipliers, enforcing the given time-variant Dirichlet-boundary conditions. 


t 
99, Q IQ, 
Os, 
Sf 4 
9% S 


Figure 1.1.: Illustration of the space-time domain for @ = 1. 


For the spatially one- and two-dimensional case, ie. « = 1 and a = 2, the space-time 
domain takes the form of a plane and cylinder, respectively (cf. Fig. 1.1 and Fig. 1.2, for a 
visualization thereof, including all relevant boundaries). 


1.1. Motivation and formulation of the problem 


> 


Figure 1.2.: Illustration of the space-time domain for @ = 2. 


Commonly initial boundary value problems are solved by applying sequential space-time 
integration methods. For that purpose, the underlying partial differential equation is 
commonly transformed into a system of ordinary differential equations by applying e.g. 
the finite element method based on a spatial discretization (see, e.g. [34, 62]). Boundary 
conditions might be taken into account either by choosing suitable test function spaces 
or by imposing feasible geometric constraints to the configuration space. The motion of 
the system is then formulated in terms of redundant coordinates on a specified constraint 
manifold. The resulting spatially discrete system can be solved by applying appropriate 


time integration schemes mostly based on finite differences”. 


In principle the same approach can be applied to the inverse dynamics problem. However, 
since the constraints do not have collocation property in the sense that the corresponding 
constraint forces do not refer to the same location as the constraints does, they need to 
be distinguished from classical geometric constraints. From a geometrical point of view, 
the constraint forces are non-(ideal) orthogonal to the constraint manifold. This specific 
type of constraints are frequently referred to as servo-, control- or program-constraints 
(see [79], [64] and [22], respectively). 


Due to this inherent non-standard construction of the constraint realization the resulting 
differential-algebraic system of equations is known to be characterized either by internal 
dynamics (see e.g. [56],[82] and [87]) or by excessively high differentiation index along 
with excessively high demands on the prescribed constraint function (see e.g. [28], [25], 
[13], [87], [5] and [63]). In the former case, the internal dynamics tend to be unstable 
depending on the spatial discretization. This hinders a numerical stable integration of the 
problem at hand. Therefore, it is inevitable to identify the unstable internal dynamics and 
carry out relevant analysis thereof (cf. e.g. [17], [85], [14]). In the latter case, appropriate 


2 Actually, instead of solving a hyperbolic problem, elliptic problems are solved successively in time 
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index reduction techniques need to be applied prior to the numerical time integration 
since a numerically stable integration of the resulting differential algebraic system of 
equations is depending significantly on the differentiation index (see e.g. [24] and [4]). 
Unfortunately it turns out, that both, the differentiation index as well as the demands on 
the smoothness of the imposed constraint increases with the refinement of the spatial 
discretization. 


The highly restrictive applicability of solving the inverse dynamics problem sequentially 
in space-time motivates this work. Our goal is to develop novel solution strategies cir- 
cumventing the issues outlined above. Therefore, analysing the governing equations in 
more detail seams to be inevitable. In particular, exposing the underlying hyperbolic 
structure will lead to a much better understanding of how information flow in physical 
time dependent problems. The causality of the solution, i.e. that the solution depends on 
the past but not on the future, will provide valuable insights into the problem and will 
turn out to be crucial for the inverse dynamics of infinite-dimensional systems. Conse- 
quently, simultaneous space-time discretization approaches may prove to be the natural 
choice for the solution of the inverse dynamics problem. This may confirm the famous 
quotation of T.J.R. Hughes and G.M. Hulbert as given by 


In fact, the more one thinks about it, the more apparent it becomes 
that the ubiquitous semidiscrete approach is conceptually confining, even 
schizophrenic. [52] 


This will pave the way of successfully developing numerical integration methods that are 
capable to solve the inverse dynamics of infinite-dimensional systems. 


Representative numerical examples including geometrically exact formualtions of slen- 
der structures undergoing large deformations, such as ropes and beams as well as gen- 
eral nonlinear continuum formulation demonstrate the capabilities of the newly devised 
space-time integration methods. Also flexible multibody systems will be addressed, un- 
derpinning the relevance of the proposed methods. 


1.2. Outline 


The rest of this work is outlined as follows. 


Chapter 2 aims to introduce various mechanical systems aligning with the initial bound- 
ary value problem postulated abstractly in Section 1.1. In particular geometrically exact 
formulations governing the motion of slender structures such as strings and beams are 
derived in Section 2.1 and Section 2.2, respectively. Section 2.3 demonstrates that the 
general theory of three-dimensional continua aligns with the postulated formulation as 
well. 


1.2. Outline 


In Chapter 3 classical solution strategies for the inverse dynamics problem that are based 
on a sequential space-time integration, are addressed. Therefore, Section 3.1 aims to give 
a brief survey of sequential space-time integration methods for the initial boundary value 
problem at hand. This will be done first in context ofthe direct dynamics problem. In par- 
ticular, the pure Neumann boundary problem will be considered first, before different pos- 
sibilities of imposing Dirichlet boundary conditions in general, are discussed afterwards. 
Based on this, the inverse dynamics problem is introduced in context of spatially finite- 
dimensional mechanical systems subjected to time-varying servo-constraints in Section 
3.2. A detailed analysis of such constrained systems aims to elaborate the fundamental 
distinctions of servo-constraints to classical contact-constraints. Consequences thereof, 
regarding the construction of numerically stable integration methods, will be addressed 
likewise. Numerous selected examples will be given. 


Due to the highly restrictive applicability of solving the inverse dynamics problem se- 
quentially in time, the governing partial differential equation is analyzed in Chapter 4 
in more detail. In particular, exposing the underlying hyperbolic structure of the gov- 
erning partial differential equations anticipates to gain more insights into the problem 
at hand. Enlighting resulting mechanisms of how information flows within continuous 
structures will pave the way to solve the inverse dynamics problem at hand numerically 
stable. Therefore, in Section 4.1 and in Section 4.2, novel numerical methods are pre- 
sented, that are based on a simultaneous space-time discretization. Several numerical 
examples are presented in Section 4.3, underpinning the relevance of the proposed meth- 
ods. 


Chapter 5 aims to apply the newly established theory of simultaneous space-time in- 
tegration to the control problem of flexible multibody systems. Therefore, a strategy to 
solve the cooperative control problem of a rigid body actuated through several flexible 
ropes is established in Section 5.1 before a rotational arm model consisting of two rigid 
and one flexible arm, is analyzed in Section 5.2. Accompanying numerical examples are 
underlining the relevance ofthe proposed methods in context of flexible mechanical sys- 
tems. 


Conclusions will be drawn in Chapter 6. 


2. Infinite-dimensional mechanical 
systems 


In Sec. 1.1, the inverse dynamics of underactuated flexible mechanical systems have been 
formulated in terms of quasi-linear hyperbolic partial differential equations subjected 
to time-varying Dirichlet boundary conditions that are enforced by unknown, spatially 
disjunct, i.e. non-collocated Neumann boundary conditions. The present chapter aims 
to derive various mechanical systems, aligning with this specific initial boundary value 
problem. Therefore, slender structures such as strings and beams as well as general three- 
dimensional continua are analyzed in Sec. 2.1, Sec. 2.2 and Sec. 2.3, respectively. 


2.1. Geometrically exact string formulation! 


This section aims to introduce a geometrically exact formulation for strings undergoing 
finite deformations, that aligns with the framework postulated in Sec. 1.1. An exhaustive 
account of the underlying geometrically exact description of strings can be found in the 
book by Antman [7]. To identify material points lying on the reference curve of the 
string, we make use of the arc-length in the reference configuration s € S = [0,L] CR 
(cf. Fig. 2.1). In this connection, we assume a stress-free reference configuration in which 
the string has length L. 


The placement of a material point s € S at any time t € T = [0, œ) is characterized by 
the deformation map r(s,t):SxT > R@, where d € {1, 2,3} is the spatial dimension. In 


the sequel it proves convenient to introduce the space-time domain Q = S x T CR. 


Balance of linear momentum gives rise to the underlying equations of motion which take 
the form of the following system of partial-differential equations 


asn(s, t) + b(s, t) = (pA)(s)a?r(s, t) (2.1) 


1 This section partly reproduces [93] 
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for all (s,t) € Q. Here, n(s,t) : Q +> Rf represents the contact force at (s,t) € Q, 
b(s,t) : Q > R? denotes the body force per unit reference length at (s,t) € Q and 
(pA)(s) : S++ R, is the mass density per unit reference length at s. In this connection, p 
stands for the density per unit reference volume and A refers to the cross-sectional area 
of the string in the reference configuration. The body force per unit reference length is 
given by 

b(s,t) = (pA)(s)g (2.2) 


where g is the vector of gravitational acceleration. 
b 

f b(s, t) ds 
N -n (a, t) a 


FR Aa 
r(s, a a A t) 


0 L 


S 


Figure 2.1.: Illustration of the geometrically exact string. 


Since the string has no bending stiffness, the contact forces in the string need to be ori- 
ented tangentially to the string. That is, 


z ðsr(s, t) 
n(s, t) = N(s, t) lern] $ (2.3) 


The tension N (s, t) at (s, t) € Q is determined by a frame-indifferent constitutive relation 
of the form 


N(s,t) = Ñ (v(s, t), s) (2.4) 
where the stretch at (s, t) € Q is defined by 
v(s, t) = |lðsr(s,t)|| . (2.5) 


The total actual length of the string follows from 


I(t) = fe t) ds. (2.6) 


Note that the stretch represents the local ratio of the deformed to the reference length of 
the string. A physically meaningful elastic constitutive law needs to satisfy the restric- 
tions N(v,s) > œ as v > œ and N(v,s) > -œ as v > 0. In addition to that, we 
assume a natural reference configuration implying N(1,s) = 0, ie. vanishing tension in 
the reference configuration. 


2.1. Geometrically exact string formulation 


Example 2.1 (Hyperelastic material model). Subsequently, we apply a hyperelastic consti- 
tutive model for a homogeneous rope based on the stored energy U(v) = A -21nv-1). 
This stored energy represents a simple model for uniaxial rubber-type material behavior. The 
tension in the string follows from N(v) = Ü’(v) leading to the constitutive relationship 


N(v) = = ot ae (2.7) 


Linearization of the last equation about the natural reference configuration yields the lin- 
earized constitutive relationship 


Min(v) =EAv-1). (2.8) 


Thus, EA can be regarded as axial stiffness in the reference configuration. 


For the inverse dynamics problem, we assume that the motion of the elastic string is 
partially specified at its boundary at s = L (Fig. 2.2). That is, the placement of the free 
end of the string at s = L is specified for all times, leading to the boundary conditions 


r(L,t)=y(t) and n(L,t)=0 VteT (2.9) 
where y : Tre R is a prescribed function of time. 


ee 


s=0 ee 


{ei} n 


Figure 2.2.: Illustration of the inverse dynamics problem of the geometrically exact string. 


The objective is to find the actuating force vector f : T +> R@ acting on the boundary at 
s = 0 causing the partially prescribed motion of the string at s = L. The corresponding 
boundary condition can be written as 


n(0,t)=f(t) VteT. (2.10) 


2. Infinite-dimensional mechanical systems 


The resulting feedforward control problem gives rise to the initial boundary value prob- 
lem postulated in Sec. 1.1. For the geometrically exact string formulation, the problem 
may be summarized as follows: Find r (s, t) and f(t) such that 


A(s)d?r(s, t) — ds (B(s, t)ðsr (s, t)) = C(s, t) Y (sH EQ, 
B(0,t)ðsr(0,t) = f(t), r(L,t) = y(t), B(Lt)ösr(Lt)=n(t) Y tET, (2.11) 
f(0)=f, and r(s,0) =ro(s), &r(s,0) = vo(s) V ses, 
where, the abbreviations A: S > R?4,B: Qh R*4 and C : Qh Rf, as given by 
A(s) = (pA)(s)Ia, B(s,t) = xen Iq and C(s,t) = b(s,t) (2.12) 
v(s, 


have been introduced. Note that Iq denotes the dxd identity matrix. The system of partial 
differential equations (2.11); results from substituting (2.3) into the equations of motion 
(2.1). Note that (2.11); can be classified as a second-order quasilinear hyperbolic system 
of partial differential equations in one dimension for r (s, t) (cf. [7]). The nonlinearity of 
(2.11); is reflected in the dependence of B(s, t) on the unknown vector r (s, t), cf. (2.12)2 
together with (2.4) and (2.5). 


Moreover, (2.11)2 corresponds to the boundary conditions (2.10) and (2.9), respectively. 
Accordingly, if the end at s = L of the string is free, we have n(t) = 0. Note, however, 
that the force vector n(t) has been introduced to take into account more general scenarios 
such as those mentioned in Remarks 2.2 and 2.3. 


Initial conditions are accounted for by (2.11)3, where fọ € R? and ro,vo € Rf are pre- 
scribed. Feedforward control problems often deal with rest-to-rest maneuvers. There, 
the solution of the corresponding equilibrium problem supplies the initial values f}, ro 
and vp = 0. In the more general case in which the initial configuration of the string 
is not at rest, the initial values for the inverse problem can be obtained by solving the 
corresponding forward dynamics problem. 


Remark 2.2 (Additional point mass). An additional point mass M might be attached to 
the boundary of the string ats = L. In this case, the last boundary condition in (2.11)2 would 
have to be replaced by 

n(t) =M (g- r(L, t)) (2.13) 
where ö?r(L,t) = D?y(t) and, as before, g denotes the vector of gravitational acceleration. 


Remark 2.3 (Multibody system). If the elastic string is a sub-system belonging to the 
inverse dynamics problem of a multibody system, the force vector n(t) in the last boundary 
condition in (2.11), is a prescribed function of time (see Chap. 5 for further details). 


Remark 2.4. In the context of infinite-dimensional models for strings, flatness-based meth- 
ods have been proposed in [55, 60, 76] to solve the inverse dynamics problem (see Def. 3.23 
and accompanying examples and references therein for more background on differentially 
flat systems). The flatness-based methods typically rely on mechanical modeling assump- 
tions such as small deformations [55, 76] or inextensibility [60]. 
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2.1. Geometrically exact string formulation 


Remark 2.5 (Linear-elastic bar). Concerning small longitudinal deformations about a 
straight reference configuration, we recover the model of a linear-elastic bar (Fig. 2.3). 


u(s,t) 
fO M 


Figure 2.3.: Illustration of the linear-elastic bar. 


Accordingly, consider the one-parameter family of configurations r° (s, t) = r° (s, t)e, where 
e is a unit vector and r° is a power series expansion of the form 


r*(s,t) = s + eu(s, t) + O (e°) . 


Here, r° = s specifies the placement of material points in the straight reference configuration 
and the derivative |, r(s,t) = u(s,t) characterizes small longitudinal displacements. 


Moreover, the resulting axial strains can be obtained via 


d 
e(st) = —| dgr*(s,t) = d,u(s,t) . 
de €=0 
Similarly, the linearized stretch yields v = 1+ d,u, so that the linearized constitutive relation 
(2.8) can be written as Nin = EAe. Eventually, linearizing the equations of motion (2.11), of 
the string about the straight reference configuration and neglecting the body forces yields 


pA u — EA u=0. (2.14) 


This partial differential equation governs the longitudinal motion of the linear-elastic bar 
and coincides with the wave equation in one dimension. The inverse dynamics problem 
(2.11) can now be restated for the linear-elastic bar: Find u : Q —> R and f : T œ> R such 
that 


u - ou = 0 YV (sHERQ, 
EA9su(0,t) = f(t), u(L,t) = y(t), du(L,t)=0 Y tET, (2.15) 
r(s,0) =s, Yır(s,0)=0 Y seS. 


Here, for the sake of clarity and without loss of generality, a free end at s = L is assumed. In 
(2.15), c denotes the constant wave propagation speed defined by 


c=4/—-. (2.16) 


There exists an analytical solution to problem (2.15) that will serve as reference for the nu- 
merical methods developed in the sequel. 
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Remark 2.6 (Analytical solution to Ex. 2.5). An analytical solution to problem (2.15) is 
based on d’Alembert’s formula (cf. e.g. [39] or [36, Ch. 16]). It can be easily verified that the 
ansatz 

u(s,t) = ®(t + cts) + Y(t — cts) (2.17) 


satisfies the partial differential equation (2.15),. Now, the first boundary condition in (2.15)2 
yields f(t) = EAd,u(0, t), where (2.17) implies 


a,u(0, t) = (D (t) — ¥’(t)) ec? . (2.18) 


Inserting (2.17) into the second boundary condition, u(L, t) = y(t), and subsequently differ- 
entiating with respect to time yields 


D (t+c'L)+¥ (t-c'L) =y (t). (2.19) 
Similarly, the third boundary condition, i.e. ðsu(L, t) = 0, gives rise to 
®'(t+c1L)-—W’(t—c!L) =0. (2.20) 


Adding (2.19) and (2.20) yields 20’(t + c7!L) = y’(t), while subtracting (2.20) from (2.19) 
yields 2¥’ (t-c"!L) = y’(t). Evaluating the last two equations att = t-c !Landt = t+c"1L, 
respectively, leads to 


£ 1 [2 = 
P) = zy (t-01), 


7 (2.21) 
Y(t) = art +L). 
Inserting (2.21) into (2.18), the first boundary condition eventually yields the result 
EA , - , - 
FO = YE- eL) - y (t+ D) , (2.22) 


providing the actuating force at the left boundary of the bar in terms of the prescribed dis- 
placement of the right boundary of the bar. 


Remark 2.7 (Analytical solution to Ex. 2.5 - alternative). Considering the one dimen- 
sional linear wave equation (2.15); governing the motion of linear elastic bars, a Laplace 
transformation eventually gives rise to an analytical solution as well. Therefore, the follow- 
ing second-order ordinary differential equation with constant coefficients 


TU(s,r) - c°? U (s, T) = 0 (2.23) 
subjected to the following boundary conditions 
EAd,U (0, t) — F(T) = 0, EAd,U(L,t) =0 and U(L,r)-T(r)=0 (2.24) 


can be obtained. In Eq. (2.24), T(r) and F(r) denote the Laplace transforms of the tra- 
jectory and actuating force, respectively. Consulting any literature on ordinary differential 
equations (e.g. [92]), the exponential ansatz 


U(s,r) = A(t) exp(As) , (2.25) 
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2.1. Geometrically exact string formulation 
solves the differential equation (2.23) (e.g. [92]). Inserting (2.25) as well as its first and second 
derivative 

d;U(s,t) = AA(r) exp(As) and 92,U(s,r) = A?A(r) exp(As) 


into (2.23) yields 
T° A(T) exp(As) — c?A* A(z) exp(As) = 0. (2.26) 


Since A, = —tce~! and Az = rc", solves (2.26), it applies that 
U(s,r) = Aı(r) exp(—rse7!) + Az (T) exp(rsc7') (2.27) 


constitutes a fundamental solution to (2.23). Taking the boundary conditions (2.24) into 
account, it follows 


U(s=L,r) = A1 (T) exp(-tLe”') + Ag(t) exp(tLe7!) =T(r), 
EAð;U (L, T) = EA (-c A, (T) exp(-tLe"') + 7 'rA2(r) exp(tLe"'))=0, 
EA9,U(0,r) = EA (-c"'rAı(r) + c"'rA;(r)) =F(r). 


The unknown functions A: (r), A2(t) and F(r) are uniquely determined as 
Aı(r) = 5 exp(rLe FL?) ; 
A(T) = 5 exp(—the P(r) ; (2.28) 
F(r) = O (exp(-tLe”') -exp(rLc"')) . 


Inserting (2.28), and (2.28), into (2.27) leads to 


U(s,r) = TO (exp(- (t(x -L))r) +exp(- (ce (L -x))r)) (2.29) 
and consequently 
a,U(s, T) = -c tA; (1) exp(-c"'rs) +c!rAz(r)exp(c'rs). (2.30) 


Transforming back into the time domain, by taking advantage of the shifting rule leads to 
the solution 


1 1 
u(s,t) = „re =c! (s — L)) + ait -cHL-S). (2.31) 
The actuating force can then be determined as (cf. (2.22)) 


ft) = = (rt - eL)-y(t+ c'L)) é (2.32) 
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2.2. Geometrically exact beam formulation 


This section aims to briefly derive the classical equations of motion of geometrically exact 
beams. Furthermore we will demonstrate, that these equations align with the framework 
proposed in Sec. 1.1. Essentially, the derivations, that are made subsequently, are based 
on [7], [88], [77] and [20]. 


Kinematics. The motion of any material points € S = [0,L] c R of a beam at any 
time t € T = [0, œ) C Ris defined by the deformation map 


r:QHR (2.33) 
determining the configuration of a curve’, along with a moving orthonormal basis 
di: Q = R? Vi € {1, 2,3} (2.34) 


indicating the average orientation of the cross-section, where d3 is normal and d; and d3 
are tangential to the cross-section. In the reference configuration, the beam is defined 


by 
rl(s,t)€a% = R(s) : u. > R? and dilist)eaa = Di(s) : AQ > R? 


for 90, = S x {0}. By introducing the proper-orthogonal tensor A € SO(3), the rotation 
of the orthonormal basis (2.34) is given by 


di = AD;. (2.35) 


At this point, it should be emphasized that planarity of the cross-section is assumed. 
Abandoning this, would require further spatial variables. 


Figure 2.4.: Illustration of the kinematics of the geometrically exact beam formulation. 


Differentiating the orthonormal basis (2.34) with respect to the spatial variable s yields 


asdi = (0,A)D; + Að; D; = (3sA)AT d; + A(O;Ao) ATAT d; . (2.36) 


2 This axis must not necessarily coincide with the line of centroids of the beam 
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2.2. Geometrically exact beam formulation 


In Eq. (2.36), the proper-orthogonal tensor Ap € SO(3) has been introduced determining 
the shape of the beam in reference configuration, ie. D; = Ape; and consequently 


sD; = Os(Age;) = (AsAo)e; + Aods€i = (AsAo)A) Dj = (sA) AGA’ d; . 


By introducing the skew-symmetric curvature matrix? 


0 —K3 K2 
Sk = (3AA =| k; 0 -K (2.37) 
—K2 Ky 0 
Eq. (2.36) can be rewritten as 
asdi = (Sk + ASLA! )d; = Syd; . (2.38) 


Note that the assumption of a straight reference configuration implies Ap = I and hence 
S? = 0. Introducing the axial vector x = x;d;, Eq. (2.38) can be defined alternatively as 


od; =K Xd; . 


Analogously, the temporal change of the moving basis can be obtained by introducing 
= widi, as 
Id; = Sodi =@WX d; : 


Herein So represents the skew-symmetric angular velocity matrix. Note, that by defini- 
tion Sow = 0. Following [7, pp. 260-261], the strain variables are defined as 


yi=oOsr-d; and xkj=K-d;, (2.39) 


where yı and yz measure shear, y3 measures dilatation, kı and kz measure flexure and x3 
measures torsion. 


Example 2.8 (Planar kinematics). Regarding planar motion, the rotation around the di- 
rector d2(s,t) = D2(s) for allt € T can be described by the following proper-orthogonal 


matrix 
cosO 0 -sin® 


A=|o ı o |. (2.40) 
sin® 0 cos® 


With (2.40) the angular velocity matrix Se = d;(A)A' simplifies to 


-sin® 0 -cos®||cos® 0 sin® 00 -1 
s,=2948| 0 0 0 0 1 0 |=90]0 0 0o|. (2.41) 
cos® 0 -sin®| I-sin® 0 cos® 10 0 


3 Orthogonality of A implies AAT = I. Differentiation of this orthogonality condition with respect to the 
spatial variable s yields ðs (A)AT + Að, (AT) = 0 which indicates the skew-symmetry of ðs (A)AT 
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Consequently the axial vector is given by œ = wzdz = W2e2. The same can be shown for the 
curvature K = K2d2 = Kze, by computing the curvature matrix in (2.37) using (2.40). 


Figure 2.5.: Planar illustration of the kinematics of the geometrically exact beam formulation. 


Equilibrium. After having addressed the kinematics, the corresponding dynamics will 
be investigated subsequently. The (material form of the) balance of linear momentum on 
an interval S D Z = [c,s] of a beam can be established as follows 


n(s,t) —n(c, t) + fae) dé; = 0;P(s,t) . (2.42) 


c 


Therein, the contact force n : Q > R?, the external load ñ : Q > R? and the linear 
momentum P : Q > R? acting on a beam segment T have been introduced. 


Since the center of gravity of each cross-section is specified as rs = r + &d, for 
E = A7! f ča dA = A™!Sg, where Sa : S > RVa € {1,2} is the first moment of area 
with respect to &, € R (cf. Fig. 2.6), the linear momentum P : Q +> R? can be obtained 
as 


P(s,t) = f ars dm 
= f "(pA)(s)a(r-+ Edy) de 
a f Gara / p(s) l ETETA 


2 f Aar rada dé vae{La 


= [ p6) an. 


Therein, the linear momentum density function p : Q +> R? as given by 


P = (pA)(s)arr + (PSa)(S)ðrda (2.43) 


16 


2.2. Geometrically exact beam formulation 


has been introduced. Defining the linear momentum density relative to r(s, t) as 
drq (s, t) = (pSa)(s)arda Va € {1,2}, (2.44) 


the derivative of the linear momentum density p(s, t) with respect to time is obtained 
as 


a:p(s.t) = (pA)(s)ar + (pSa)(s)ada = (pA)(s)aer+ aq Va e {1,2}. (245) 


Obviously, the relative linear momentum density function q : Q > R, defined in (2.44), 
vanishes by choosing r(s,t) accurately, ie. r(s,t) = rs(s,t), such that Sg = 0 for 
a € {1,2}. Analogously, the material form of the balance of angular momentum for a 
segment Z can be established as 


m(s,t) — m(c, t) + (r(s,t) x n(s,t)) - (r(c, t) x n(c, t)) 


s s 2.46 
i J ELORE 1 ater. © 


Therein, the contact torque m : Q + R°, the external applied torque m : Q +> R? and 
the angular momentum (relative to an inertial frame) L : Q +> R? of a beam segment T 
have been introduced. Taking advantage of the definition of the second moment of area 
Iag = F Čačp dA and assuming rp = r+&.d, (cf. Fig. 2.6), the angular momentum L(s, t) 
can be obtained as 


L(s,t) = / rp X o:(rp) dm 
— fe + čada) X (Orr + Epd:dp) dm 
= e f T rx ardas f &p(r X ðıdp) + alda X Hr) dA 
ö (2.47) 
+ [ tačplda x dp) dAd, 
= / pA(r x Or) + PSa(r X da + da X Ir) + plap(da X Adz) dé 


= [un dé; Va, p € {1,2}, 


where the angular momentum density function I : Q +> R? as given by 
I(s,t) = pA(r X Gr) + pSalr X Ada +da X Ar) + pla(da X Ada) Væ, fp € {1,2} (2.48) 
has been introduced. By defining the angular momentum relative to r (s, t) as 


h(s, t) = plap(da X dg) Va, B € {1,2} (2.49) 


2. Infinite-dimensional mechanical systems 


the time derivative of the angular momentum density function can be established as 


ıl (s, t) = pA(r x 8r) + pSalda X Or) + pSa(r X da) + Plap(da X adp) (2.50) 
= pA(r x dr) +rxdqtqxart+oh. l 


Figure 2.6.: Illustration of the geometry of the cross-section of the geometrically exact beam. 


Taking the derivative of (2.42) and (2.46) with respect to the spatial variable s € S, the 
balance equations are represented in local form as given by 


ðsn +A = dp (2.51) 
a.m + (r X dn) + (ar Xn) + (r Xn) +m = Al. 


Inserting Eq. (2.45) into Eq. (2.51)ı, the following relation can be established 
rX op = pA(r X der) +r X 0eq=rxXdntrxn. 


By introducing al = qx ö?r+9;h, the balance of angular momentum in local form (2.51)2 
might be reformulated as 


am +oasrXn+m = ol, (2.52) 


Equation (2.51); and (2.52) are the equations of motion for (the classical form of Cosserat) 
rods ([30], [7]). These equations are also frequently and synonymously referred to as the 
geometrically exact ([20]) or the Simo-Reissner beam formulation ([88],[89], [77] and [78]). 
It might be convenient to write them compactly as 


fl,- lal,- 
I st m „Ss 
Example 2.9 (Planar equilibrium). For the plane case, the angular momentum relative to 
r(s,t), defined in (2.49), and its time derivative reduces to 


0 
ðr xn 


+ 4 l (2.53) 


m 


h(s, t) = plıı(dı x dıdı) = w2e2 and ðth(s, t) = (3rw2z)ez, 
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2.2. Geometrically exact beam formulation 


respectively, due to d;dı = dz = 0 and 
dı x ôrdi = Ae; X S„Aeı = W2€2. (2.54) 


Therein use of the kinematical relations, presented in Ex. 2.8 for the plane case, has been 
made. In Fig. 2.7, the geometry of a cross-section for the planar case is depicted. 


Figure 2.7.: Planar illustration of the geometry of the cross-section of the geometrically exact beam. 


Furthermore, the contact force and contact moment are assumed as n = (nı, 0, n3)T and 
m = (0, mz, 0)!, respectively (cf. Fig. 2.8). The same holds for the external applied force and 
moment ñ and m, respectively. 


Figure 2.8.: Planar illustration of the equillibrium of the geometrically exact beam formulation. 


Subsequently, it will be demonstrated, that the classical equations of motion for Cosserat 
rods, introduced in (2.51); and (2.52), align with the mathematical framework postulated 
in Sec. 1.1. For this, the contact force and moment can be written alternatively as 


n= N;d; = N;ANe; and m= M;d; = M;ANe; š 


Additionally we need to resolve the geometric constraints, that are imposed on the com- 
ponents of the directors accounting for the orthonormality. One possible way of doing 
this is by introducing Euler angles © : Q +> R?. The Euler angles form a local basis in 
SO(3) such that the overall rotation A(®) € SO(3) of a moving frame {d;} relative to an 
inertial frame is defined by three rotations that are carried out successively: (i) rotate by 
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©, around es, (ii) rotate by ©; around h and (iii) rotate by ©; around d3. See Fig. 2.9 for 
a graphical illustration of these rotations. 


Following [7, pp. 300-301], the strain variables are then given as 


kk = (RA,O), and wr = (Rd;O), (2.55) 
where 
sn®; —cosQ@3;sin@, 0 
R= | cos ©; sin ©; cos®; 0 
0 cos ©; 1 


Note that the introduction of the Euler angles may lead to polar singularities. For more 
details on Euler angles see [11, pp. 148-150], [70, pp. 493-495], [44, pp. 150-154], [59] and 
[69]. 


Figure 2.9.: Rotation of a moving frame {d;} relative to a fixed frame {e;} by using the Euler angles. 


Focusing on hyperelastic materials, the constitutive relations are governed by the stored 
energy function Y = ¥(y, x). We assume that 


N; = 9 ¥(y, K) = N;(y, K) = Fik(y, K)Yk > 


ý (2.56) 
Mi = I, ¥ (y, K) = M;(y, K) = Gik(y, K)Kk 


holds. Note that the fundamental conditions regarding the limiting deformation cases 
have to be fulfilled. Hence, for ya — {+00} the contact force Ng should tend to +00 
and the contact force N3 should tend to +oo for y — {00,0}. Since the principle of 
impenetrability of matter must not be violated, the contact moments M; should tend to +00 
as the curvature x; tends to an upper or lower bound, where an intersection of adjacent 
cross-sections is imminent (cf. Fig. 2.10). 
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Figure 2.10.: Limiting case defined by the principle of impenetrability of matter: Imminent intersection of adja- 
cent cross-sections. 


Example 2.10 (Planar constitutive model). Considering planar motion, the following strain 
energy function of the form 


1 EA 
Y= 7 EIk? + GAy? + ae — 2Iny3 — 1) (2.57) 


might be chosen, leading to the constitutive relations for the contact forces N; and N; as well 
as for the contact moment M} as given by 


EA 
N; = ie -y;'), Nı=GAyı and M,=EI x2. (2.58) 


Note that the following fundamental conditions need to be ensured 


N3(ys) > +00 for y > ey 
Ni(y1) >+% for yı > +% (2.59) 
Me (K2) +0 for mo Be {k2 : y3 > V(k2,s)} , 


where condition (2.59); ensures a non-intersection of neighboring cross-sections (cf. [7]). 


By taking the constitutive relation (2.56) into account and recalling (2.39), the contact 
force and torque might be reformulated as 


n =A- Fiklek 8e;)A! - dar = (AFTA’) - or (2.60) 


and 

m = AGix(ex Q e;)-R'0,@ = (AG'R’)-0,0, (2.61) 
respectively. In Eq. (2.60) and Eq. (2.61), F = Fiz(y)(ei ® ex) G = Gix(y)(e; ® ex) has 
been introduced and use of the identity Ab @ Cd = [C(d ® b)AT|" has been made’. 
Recalling that d;d. = w X da = Soda, the second time derivative of the directors d can 
be established as 


Bda = 0;(Soda) = :(Sw)dq + Sda = (Hw X da) +S? da = SZ da —dg X Hw . (2.62) 


4 This property follows directly from the definition of the dyadic product (a - b)c = (c @ a) - b (cf. [50, p.10]) 
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This leads to the time derivative of the linear momentum relative to r(s, t) as given by 
A:q(s,t) = (pSa)(s) (-da X dro + Soda) 
= 5,94 - (q X 40) (2.63) 
= 5,9 - Sq: (R2) . 


Note that Se and Sg, introduced in (2.62) and (2.63), respectively, denotes the skew- 
symmetric matrix 


0 -()3 O P 
Soy =] (+) > Sr]: 
-()2 Ch 0 


Taking advantage of this relations, the balance of linear momentum, Eq. (2.51); can be 
reformulated as 


(pA)(s) ar — SqRA2O - ð, (AF™ATa,r) =n — S2 q +$q0,Rd,0 . (2.64) 


Furthermore, by consulting (2.49) for a definition of the angular momentum relative to 
r(s, t), its time derivative might be rewritten as 


Ah(s, t) = (plap)(s)dı (da X dp) 
= (plap)(s) [da x (31% x dg) + da X (@ X (w x dp)) + (@ X da) x (w x dp)] 
= (plap)(s)[da X (O(R9,O) x dg) + da X (@ x (w x dg)) 
+ (w X da) X (w x dg)], 


such that the balance of angular momentum (2.52) can be reformulated as oe 
Saar - (pJ)(s)aO — ð, (ATR ; 2,0) -arxn+m-(pk)(s), (2.66) 
where the functions J : Q + R33 and k : Q + R? have been introduced as? 
JO) = (lep)(s)Sa,Say (2.67) 
and 
k(s) = (lap) (s) (da X S? dg + @ - (da X dg)@) + dg X (4, RO x dg)) . (2.68) 
In the context of the inverse dynamics problem, the actuating force and torque 
fm: dQ Rt (2.69) 
applied to 90, = {0} x T, causing a partly specified motion 
XI(s2)ea0, = Y(t) : 90, > RI. (2.70) 


5 Note, that use of the (a x b) x (a x c) = a- (b x c)a has been made. 
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on dQ, = {L} x T are searched for. (cf. Fig. 2.11 for a planar illustration of the inverse 
dynamics problem) 


Figure 2.11.: Planar illustration of the inverse dynamics problem: Find the force f(t) and moment m(t) acting 
at s=0 such that the imposed time-varying constraint on d;(t) and O(t) ats = L is met. 


Regarding multibody systems, Neumann boundary conditions in the form of a system of 
ordinary differential equations 


Bdsx|(s1)ca0, = (x, tx, x, t) (2.71) 
may be taken into account on dQ, = {L} x T. By specifying initial conditions on 
OQ) =S X {0} 

x| (5,1) €aQ) = Xo IX|(s1)eoo, = Vo (2.72) 


and introducing the coefficients 


ñ — S},q + 599, R9,O 
(dsr Xn) +m — pk 

(2.73) 
the problem aligns with the general formulation of the inverse dynamics of flexible me- 
chanical systems proposed in Sec. 1.1. 


As a || Ge ea 0 


Sq -pJ 0 nor”) on al 


Example 2.11 (Planar problem). Subsequently, the special case of planar motion by ad- 
ditionally assuming straight reference configurations as well as r = rs is considered. Then 
strain energy functions of the form (2.57) can be assumed (cf. Ex. 2.10) such that the functions 
F(y) and G(x), introduced in (2.60) and (2.61), respectively, appear in diagonal form, i.e. 
F = diag(GA, 0, zA (1-y3°)) and G = diag(0, EI, 0) . Consequently, the diagonal block 
matrices AF’ A’ and AGTRT in (2.73), are obtained as 


cos? © 0 cos®sin® 
AFTA! =GA 0 0 0 
cos®sin® 0 sin? © 
EA sin? © 0 -cos®sin® 
+— (1-y°) 0 0 0 
2 > 
-cos®sin® 0 cos’ © 
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and 
0 0 0 
AG'R=|0 EI ol. 
0 0 0 


Note that the Euler angles ©, and ©; vanishes in the plane case. A has been taken from 
Ex. 2.8 such that the coefficient B in (2.73) can be identified as 


cos? © cos®sin® 0 
EA _2 À 2 
B=—(1-v “)|cosOsinO sin’ © 0 
0 0 0 
GAsin?® -GAcos@sin20 0 
+ ]—-GAcos@ sin 20 GA cos? © 0 
0 0 EI 


Since the functions J and k, introduced in (2.67) and (2.68), respectively, reduce to 


0 0 . 
J = phis, = pli $ -(d, + d*,) ==1 $ and k= 0, 
: 0 0 


the two remaining coefficients A and C in (2.73) can be obtained as 


pA 0 0 ny 
A=|0 pA 0 and C= ñs (2.74) 
0 0 pI Mz + (ðsr3nı — ðsrınz) 
Therein, use of 
Osrı nı 0 
ðrxn=| 0 Ix|0| = |dsranı — Osrınz 
Osr3 n3 0 


has been made. 


Remark 2.12. (Componential representation) The classical equations of geometrically ex- 
act beams have been introduced in (2.51); and (2.52) in vectorial representation. However, in 
certain cases, a componential representation of these equations might be more convenient’. 


The contact forces and torques might then be established as 
n= Nid; and m= Mid; 3 (2.75) 
respectively. By differentiating Eq. (2.75), with respect to the spatial variable s 


In = (0,N;)dj + N;ösd; = (0;.Nj)di + N;(k x d;) y 


é The terminology vectorial and componential representation for distinguishing the chosen basis can be found 
in [7]. 
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the balance of linear momentum can be stated as 
osNiöi; + Ni(« x d;)-dj+n-dj=aAp-d;. (2.76) 


Taking advantage of the property of the triple product (x x d;)-d; = « - (di X dj), Eq. (2.76) 
might be reshaped as 
O;Ni6;; + NiQij; +n- dj =p: dj (2.77) 


by introducing 
Qij ZK: (d; x dj) $ 
Differentiating Eq. (2.75)2 with respect to the spatial variable s 
om = (0;M;)d; + M;(K x d;) 


and taking into account, that” 
dsr Xn = (S’N)d; 


holds, the balance of angular momentum might be established in componential form as 


sMiðij + MiQij + N;S7, +m: dj = Al : dj. (2.78) 


Considering that ñ - d; = (Aei) = (ATñ)e;, the equations (2.77) and (2.78) constitute the 
equations of motion of geometrically exact beams in componential form, as given by 


Ni Qij 0 Ni n:d; = ap- dj 
o : S A H t a ~ E i = : (2.79) 


Example 2.13 (Componential representation). Regarding planar motion, the functions Q 
and S? , introduced in (2.77) and (2.78), reduce to 


0 0 ke 0 -y 0 
Q = 0 0 0 and S = Y3 0 -yı 
—K2 0 0 0 yı 0 


Ni 0 K2 0 N] -dı ap: dı 
N3| +|-K2 0 O}|N3|+]a-d3| = |a%p-ds| . (2.80) 
M2 yo -yı 0j|M md; O,l - dz 


Using the hyperelastic constitutive equation, proposed in Ex. 2.10, a relation for the forces 
and torques to the strain variables can be established as 


= GAcos® GAsin® e 
Menr PAPA u -4(1 = y”) sin © Aa — ¥3°) cos® for ie {1,3} 


7 This can be shown by simply computing dsr xn = ygdx X Nidi = (yidit--+)x(Nidit---+) = (yXxN)d; 
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and 
M;i = Gijx;j = EId,O for r= 2. 


This leads to the characterizing coefficients 


cosQ sinO 0 0 kz O||N, -dı 
A=pA|-sinO cos® 0 and C=|-k. 0 O;||N3/+]n-d; 
0 0 pli Y3 = 0 Mp m: d 
as well as EA 3 a EN 
Z—(1-v”)cos® Z(1-v“)sin® 0 
B= -GAsin® GAcos® 0 
0 0 EI 


Hence, the problem at hand again coincides with the quasi-linear hyperbolic partial differ- 
ential equation, proposed in Sec. 1.1. 


Remark 2.14 (From Simo-Reissner via Timoshenko-Ehrenfest to Bernoulli-Euler). Intro- 
ducing the function u: Q > R denoting the displacement relative to the reference configu- 
ration R(s), the spatial configuration of the beam can be obtained as 


r(s,t) = R(s) + u(s,t). (2.81) 


Assuming straight reference configurations, i.e. R(s) = s e3, the derivative of (2.81) with 
respect to the spatial variable s can be determined to be ðsr = e3 + Osu. Then the strain 
variables yı and yı can be established as 


ys = ðsr - dz = (Osu + ez) - dz = (Osu + ez) - (cos Oe; + sin Ge) 


(2.82) 
= (0,u3 + 1) cos O + d,u, sin O = I; +1 
and 
yı = 9r - dı = (Osu + €3) - dı = (0, + ez) - (- sin Oe; + cos Oe, ) (2.83) 
= —(0,u3 + 1) sin © + ðsuı cosO =T}, i 
respectively (cf. [88, Eq. (5.2)]). Assuming inextensibility of the beam, i.e. 
y3 = ðsuz cos O + cos © + du, sinO = 1, (2.84) 
implies a relation for the spatial derivative of the displacement u3, given by 
1 i 
ðsuz = —— (1- cos O - d,u; sin O) . (2.85) 
cos © 
Inserting (2.85) into (2.83), the strain variable yı is obtained as 
sin® sin? © 1 > 
y=- + osu1 +cos®| = (Hu — sin O) . (2.86) 
cos © cos © cos © 
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Inextensibility of the beam implies that the force N3 is not defined constitutively anymore. 
Instead N; has to be considered as a Lagrange multiplier enforcing condition (2.84). For small 
deformations it is feasible to assume N; = 0. Then Eq. (2.80) immediately reduces to 


O;Ny +ñ = OP and sM + Ny + Mz = ab 3 (2.87) 


The system of partial differential equations (2.87) can be traced back to the pioneering work 
of S.P. Timoshenkho and P. Ehrenfest [95, 96]. 


Figure 2.12.: Illustration of the strain components of the geometrically exact beam formulation and its physical 
interpretation (cf. [91]). 


A further restriction, pertaining to the rigidity of shear deformation, i.e. yı = 0 and simul- 
taneously neglecting rotatory inertia effects, leads to the classical equation of motion 


05 (3s M2 + m2) +ñ = OP $ (2.88) 


Eq. (2.88) has been developed originally by L. Euler, Jas. Bernoulli and D. Bernoulli, by 
additionally assuming linear constitutive relations of the form 


M = Elk = Eld,@ = Ela? . (2.89) 


2.3. General nonlinear continuum formulation 


After having considered slender structures such as strings and beams in Sec. 2.1 and 
Sec. 2.2, respectively, the elastodynamics of general three-dimensional continua will be 
considered in this section. Therefore, we aim to elaborate briefly the cornerstones of the 
underlying theory including large deformations. Analogously to the theory of slender 
structures, it might be convenient to start with a recapitulation of the classical concepts of 
the kinematics of continua, before fundamental balance equations, governing the motion 
of arbitrary deformable bodies B, can be introduced. 
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Kinematics. General a-dimensional continua might be perceived as arbitrary shaped 
and deformable bodies occupying regions in a-dimensional space, i.e. B c R“. We will 
focus on formulations based on material points s € B at time t € T, i.e. material coordi- 
nates (s,t) € Q are chosen as independent variables. Such formulations are frequently 


referred to as material formulations, or more historically Lagrangian formulations?. 


The position of every material point s € B for every time t € T is then defined by the 
injective? map 
(s,t) E Q} r(s,t) ER”. (2.90) 


Similar to the one-dimensional case of slender structures such as strings and beams, it 
is inevitable to find feasible measurements in order to quantify the deformation of the 
body. Therefore, we will investigate the change of a curve so(&) : Ee [a,b] > RY, 
parameterized by the arc-length parameter £ (see Fig. 2.13). 


der (So, t) 


r (So, t) 


Figure 2.13.: Kinematics of a general non-linear continuum formulation. 


By defining the deformation gradient F = dsr(so,t), the vector field tangential to the 
curve So is obtained as 


der (So, t) = Asr(So(E), t) - Deso = F(So, t) - Deso . 


In order to describe the change of the shape of the curve s locally, we can compare the 
length of the curve so after some deformation 


b b b 
i= [ örr (so, ©) || a= f llösr (so, t) - Des| a= f ||F(so,t)-Des|| aë (2-91) 


to its length in reference configuration 


b 
b= f [Deso de 


o0 


Alternatively, the problem can be formulated in terms of the positions r (q (y, t), t) of material points cover- 
ing positions y at time ¢, i.e. spatial coordinates (y, t) are chosen as independent variables. Therefore, this 
formulation is frequently referred to as spatial or historically Eulerian formulation. 

Injectivity of r (s, t) needs to be supposed, due to the principle of impenetrability of matter, i.e. two distinct 
material points cannot simultaneously occupy the same position in space. 


Noy 
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and analyze the limit case of its ratio, i.e. 


ba 


l 1 
lim cs ||F(so,t) - des] = [(F - Deso) - (F - Deso) |? 
= [Deso  F'F - Des]? eee) 

u [Des(a, t)-C- Dys(a)]? ; 


Thus, a fiber in point a, that is oriented tangentially to the curve s changes its length 
relative to its reference length according to a scalar parameter 


1 
v= [Des(a, t)-C- Des(a)| Zur 

The curve s is stretched or compressed locally in point a as v > 1 and v < 1, respectively. 
For v = 1, the length of the curve does not change (cf. Sec. 2.1, for conspicuous analogy). 
In (2.92) the symmetric and positive definite tensor C = FTF has been introduced. This 
tensor is frequently referred to as right Cauchy-Green deformation tensor. 

Based on the right Cauchy-Green deformation tensor C, a vast variety of different strain 
measurements can be established. One prominent member of this family of strain mea- 
surements is the material strain tensor 


1 
E=>(C-D 


also referred to as Green-Lagrange strain tensor. Its virtue lies in the fact that it vanishes 
in the reference configuration. 


Equilibrium. In analogy to Sec. 2.1 and Sec. 2.2, the two fundamental balance equa- 
tions regarding linear and angular momentum are introduced subsequently for the gen- 
eral continuum formulation. 

Balance of linear momentum. According to Newton’s second axiom, the resultant force 


f : Qt R” as given by 
f= | afisa) 
B 


acting on a body B with total mass 


ms) = f dm(s, t) 


is equal to the temporal change of the linear momentum of the body as given by 


= farts t) dm(s,t) — f(t) =0. (2.93) 
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Furthermore, we assume, that a decomposition of the resultant force f into body forces 
b : Q + R“ per unit of reference volume and surface traction t : Q + R” per unit of 
reference area OB, i.e. 


f= [po t) do(s, p+ fe t(s,t) da(s, t) (2.94) 


can be applied. Eq. (2.93) along with Eq. (2.94) represent the balance equation of linear 
momentum in its integral form. 


Local form of the balance of linear momentum. With the help of Cauchy’s theorem the bal- 
ance equations can be established in local form. This fundamental stress theorem, postu- 
lated by Cauchy, states that there exists locally a second-order tensor field P : Q =œ R** 
such that 

t=P:v (2.95) 


holds. Therein, v denotes the unit normal vector on the surface, on which the traction 
t acts. To prove Eq. (2.95), the following geometrical considerations are made: Consider 
a tetrahedron, that is bounded by four planes with normal vectors ex for k € {1,2,3} 
and v such that v-e, > 0 Y k € {1,2,3}. The size of the tetrahedron is defined by hv, 
where h € R*. The area of the plane II that is orthogonal to v can be obtained as ch’. 
Accordingly the planes Ix, that are orthogonal to ex have the area ch’v - ex. This is the 
projection of the area of plane II onto plane IIx. Cauchy’s theorem can then basically be 
proven by establishing the balance of linear momentum for this specific tetrahedron as 
given by 

3 


io t(s,t,v) da(s) + 2 


k=1 


t(s,t,—e;,) da(s) =f pdir(s,t) - f(s, t) do(s). 
OB, (h) B(h) 


(2.96) 
Dividing (2.96) by ch? and evaluating the limit for h > 0 gives for every material point 
rise to 


3 
t (So, t,v) =- X Elso, t, —ek) 8er) V 
k=1 
and Cauchy’s theorem (2.95) is proven to be right, by introducing the second-order tensor 


field ; 
P(s,t) = - $ (ts, t,—ek) Ber). 
k=1 


This second-order tensor field is classically referred to as first Piola-Kirchhoff stress tensor 
field. Due to its lack of symmetry it sometimes proves convenient to introduce the second 
Piola Kirchhoff stress tensor as S = FP or the Cauchy stress tensor as ø = (det F)'PFT. 


Inserting Cauchy’s theorem (2.95) along with (2.94) into (2.93) and applying the divergence 
theorem, a local form of the balance of linear momentum is given by 


div PT +b = par. (2.97) 
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Balance of angular momentum. Besides the balance of linear momentum, the balance 
of angular momentum represents the second fundamental principle for elastodynamics, 
ie. 


= tl r(Q) x ar(Q) dm(s) - m(t,0)=0. (2.98) 
B 


Therein, m(t,0) denotes the resultant torque on the body B about the origin 0 at time t 
as given by 


m(t,0) = f rx af(a)= f r(a) xbdo(s) + f tO) dats) (2.99) 


Analogously to the balance of linear momentum, the balance of angular momentum can 
be described in its local form too. This implies 


F-P'=P.F'. (2.100) 


This might be proven by applying Cauchy’s theorem to the balance equations (2.98) and 
(2.99) straightforward obtaining 


f rx(P-v) dar f (rx f- prx ir) d=0. (2.101) 
oB B 


Multiplying (2.101) with the cross product of two arbitrary but constant vectors a € R” 
and b € R®, Eq. (2.101) reduces to 


(b8a-a®b): P- FT + (div (PT) + f - pdr) @r du=0. (2.102) 
B 
Keeping in mind that B is arbitrary and using the local form of the balance of linear 
momentum (2.97), the following relation 
(b®a-a®b):P-FT 


is obtained and the symmetry postulated in (2.100) is proven to be true. 


The partial differential equation (2.97) along with (2.100) constitute the material form 
of the classical equations of motion for general continua and aligns with the problem 
formulation postulated in Sec. 1.1 for 


A= pl, Bösrr=P and C=b 


With (2.97) and (2.100), six equations could have been derived for the twelve components 
of the unknown functions r and P. In order to determine the full system of equations, 
additionally constitutive equations need to be developed. In fact, these equations estab- 
lish a relation between the two functions P and r induced by material behavior, E.g. in 
context of elasticity the aim is to find constitutive equations of the form 


P(s,t) = P(dgr,s) . (2.103) 
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Remark 2.15 (Geometrically exact beam formulation in the context of the general non- 
linear continuum formulation). The geometrically exact beam formulation analysed in 
Sec. 2.2 can be perceived as a one-dimensional continuum. In this remark we aim to work 
out the links to the general three-dimensional nonlinear continuum formulation (cf. [88, 
App. A]). For that purpose we restrict ourself to straight reference configurations defined by 


material coordinates s = [á & gf. In Fig. 2.14 an illustration of the beam is depicted. 


n(&3,t) 


Figure 2.14.: Material (left) and spatial (right) configuration for the geometrically exact beam formulation. 


The resultant contact force per unit reference length, acting on cross-section A, in spatial 
configuration can be established as 


nn = f t(s,t) dA . (2.104) 


Cauchy’s theorem (2.95), states that the stress vector t : Q +> R? per unit reference area Ay 
acting on the cross section A, is determined by the second-order tensor field P : Q > R?’ 
and the normal vector v of the contact plane, i.e. t(s,t) = P(s,t) - v. Since, in context of 
the beam theory, the orientation of the contact surface is obviously always defined by the 
normal vector e3, the following applies 


P -e3 = (t(s,t, ex) Q ex) - e3 = (ex - e3)t(s, t, -ex) = t(s, t, —e3) . 


Therein, the definition of the first Piola-Kirchoff stress tensor field, i.e. P = t(s,t, ex) ® ex 
where k € {1,2,3}, as well as the orthogonality condition e; - e; = ôi; have been exploited. 
Eq. (2.104) can then be reformulated as 


n(é3, t) af P -e3 dAo =f t(s, t, —e3) dAo. (2.105) 


Taking into account that the divergence of the first Piola-Kirchoff stress tensor field is ob- 


tained as 
divs P = dz, Pikei ® ex = dz, t(s, t, —ex) (2.106) 
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as well as making use of the local form of the balance of linear momentum (2.97), the deriva- 
tive of the resultant contact force (2.104) with respect to & can be established as 


ðgn = / dz, (t(s, t, -e3)) dAo 
Ao 
= Ni div P — (dz,t(s, t,—eı) + dz t(s, t, -e;)) dAo 
Ao 
= i per -b- (dz,t(s, t, —e1) + dz,t(s,t, -e;)) dA, . 
Ao 
Taking advantage of the divergence theorem 
[rs t,e;) a t(s, t, —e;) -v; dT (2.107) 
A ðA 
and introducing, in accordance to (2.94), the resultant force 
f(s,t) = t(s, t, -e;) - vi ars f dA, (2.108) 
aA A 


the balance of linear momentum for the geometrically ecact beam is obtained as 


dgn = (pAo)(s)d%r(s,t) - f(s, t) . 


Compare with Eq. (2.51);. Analogously, the balance of angular momentum of the geometri- 
cally exact beam formulation can be derived from the three-dimensional theory (cf. [88]). 
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3. Sequential space-time integration 


Initial boundary value problems that occur in non-linear structural dynamics are com- 
monly solved by applying a sequential space-time discretization, cf. e.g. [34, 62]. The so- 
called method of lines is therefore typically based on a discretization in space by means of 
finite elements, followed by an appropriate discretization in time mostly based on finite 
differences. In Sec. 3.1, a brief survey of such sequential space-time integration methods 
for the initial boundary value problem at hand will be given. This will be done first in the 
context of the direct dynamics problem. For that purpose, the pure Neumann boundary 
problem is introduced before different possibilities of imposing Dirichlet boundary con- 
ditions in general, will be discussed afterwards. Based on this, in Sec. 3.2, the inverse dy- 
namics problem will be introduced in the context of spatially discrete mechanical systems 
subjected to time-varying servo-constraints. A detailed analysis of such constrained prob- 
lems aims to elaborate the fundamental distinctions between servo-constraints and classi- 
cal contact-constraints. Consequences thereof, regarding the construction of numerically 
stable integration methods, will be addressed likewise before numerous selected exam- 
ples will be given. This section partly reproduces [93]. 


3.1. Direct dynamics problem 


Considering for the moment the formulation of the classical direct, pure Neumann prob- 
lem! associated with the initial boundary value problem (1.1), then the aim is to find the 
overall motion of the system affected by external loads, neglecting any Dirichlet boundary 
conditions. 


Multiplying Eq. (1.1) by sufficiently smooth test functions w : S > R@ and subsequently 
integrating over the spatial domain S yields 


[ro - A(s) - 0?x(s, t) ds - [ro - divs (B(s, t) - 9x(s,t)) ds = [ro -C(s,t) ds. 
S S S 6 


1 The pure Neumann problem is obtained by removing the displacement boundary conditions. 
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Applying integration by parts to the second integral on the left-hand side and taking 
into account the Neumann boundary conditions leads to an equivalent weak form of the 
boundary value problem at hand 


[ro - A(s) - 0?x(s,t) ds + fowo - B(s, t) - 9,x(s,t) ds 
S S 


= [w(s) CED ds- FO wis Dlenes ËP 
+n(t)- w(s,t)|(steas, - 


This weak form can be discretized in space by applying suitable finite element approxi- 
mations to the vector-valued test function w(s) and the trial function x(s, t). We confine 
our attention to piecewise linear approximations based on Lagrangian shape functions 
(cf. [51]). The corresponding interpolations read 


w"(s) = Las) and x"(s,t) = 3 L(s)q;(t) . (3.3) 
i=1 j=l 


where L;(s) are a-linear Lagrangian shape functions and ns € N denote the number of 
nodes. The associated nodal values w; : T > R@ of the test functions are arbitrary, while 
qj: TP R? denotes the nodal position vectors at time t € T. That is, q;(t) = x"(s;, t). 


Inserting the finite element approximations into weak form (3.2) yields the semi-discrete 
equations of motion 
MD;q=F(q,t), (3.4) 


where 


F(q,t) = f(t) - FA - Br f(t) + Bon). (3.5) 


Here, q : T +> R@”s is the nodal configuration vector that contains the nodal position 
vectors at time t € T of the discrete problem at hand, i.e. 


q(t) 
gs 2% (3.6) 
qpa (t) 


Furthermore, the nodal contributions to the mass matrix M, the internal force vector 
f™ (q) and the external force vector f*'(t) are given by, respectively, 


Mij(t) = 1a | AOL ds, 

S 

f™(q.t) = f D,Li(s)B(s, t)ðsr”(s, t) ds, (3.7) 
S 

FAO = f LOCEN ds, 
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where Iq is the d x d identity matrix. The matrix By in (3.4) is of Boolean type and 
essentially links the force f(t) to the nodes lying on 907, accounting for the virtual 
work contribution dq - B; f(t) = wı : f(t) emerging from the right-hand side of weak 
form (3.2). Here, ôq contains the nodal values w; in analogy to (3.6). 


The semi-discrete equations of motion (3.4) constitute a system of non-linear ordinary 
differential equations of second-order. In the standard forward dynamics problem, the 
forcing terms f,,,(t), f(t) and n(t) are prescribed functions of time and (3.4) can be 
solved by applying common time-stepping schemes. 


Example 3.1 (Linear elastic bar). In the case of the linear-elastic bar (Ex. 2.5), the semi- 
discrete equations of motion (3.4) boil down to a system of ordinary differential equations 
with constant coefficients. Confining to piecewise linear approximations, the nodal contribu- 
tions to Eq. (3.5) can be evaluated to be 


F(t) = Kiju;(t) — 0 — ôn f (t) + ôi prn (t) (3.8) 


where for each element 


(e) _ EA 
Kr = T (51; — dr, — &i-1,) 
= 


holds. By assuming lumped masses, the nodal contributions to the mass matrix read 


(e) (e) 
OL = PEN 

My? = pA=— = — ô; . 
For the more general case of as deriving a consistent mass matrix, the contributions to the 
mass matrix are obtained as 


@_ pn ahs” 
M;; =pA A (2ôi; + Ôi+1,j + 5;-1,;) 


again by confining to piecewise linear Lagrangian shape functions. 


Remark 3.2 (Additional point mass). If an additional point mass M is attached to the right 
end of the string (Rem. 2.2), n(t) = M (g - ö?r(L,t)) has to be taken into account in weak 
form (3.2). Correspondingly, in the semi-discrete formulation, the following entries of mass 
matrix (3.7), and external load vector (3.7)3 need to be modified according to 
Mos1,p+1 — Mp+1,p+1 + MIa , 
fou Son + Mg - 


Dirichlet Boundary conditions. Boundary conditions pertaining to the configura- 
tion space Q = R@”S of the underlying mechanical system, i.e. 


x (seo = X(S 0), (3.9) 
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are called Dirichlet boundary conditions. Classically, such conditions are taken into ac- 
count by choosing suitable spaces for the test functions, i.e. 


W ={w:9 m Rf |w=00n3SpxT}. (3.10) 


The system at hand can then be formulated in a generalized, minimal configuration space, 
ie. q : T |œ R@{Rs“nD), representing the nodal configuration vector that contains the 
nodal position vectors, except for the nodes lying on the Dirichlet boundary Sn, at time 
t € T of the discrete problem at hand. The solution is then governed by ordinary differ- 
ential equations of the form 

M:D’g+F=0. (3.11) 
Alternatively, boundary conditions of the form (3.9) might be ensured by imposing fea- 
sible geometric constraints to the configuration space Q. The underlying system can 
then be formulated in terms of redundant coordinates q on a specified constraint mani- 


fold C = {q : g(q) = 0}, such that the motion is governed by a semi-explicit system of 
differential-algebraic equations” given by 


M D?q + F(Drq,q,t) +G' A(t) =0, 


g(t,q) = Gq- y(t) =0, (3.12) 


where A : T > R? are the Lagrange multipliers enforcing the given constraints g : 
TXQ => RD, Furthermore, the constraint jacobian G = Ogg (t, q) has been introduced 
in Eq. (3.12). For more details on constrained mechanical systems, see Rem. 3.12, Ex. 3.13 
and references therein. By introducing the velocity v = D:q, the system (3.12) might be 
reformulated as a first-order system 


Div = F,(0,q,A) = M! (f, q) - G71), 
Diq = F2(v,q) =0, (3.13) 
0 = F3(q) = 9(9). 
with non-singular 
OqF 3+ oF ` d Fı = -GM 'G" . (3.14) 
Then Eq. (3.13) can be identified as a system of differential-algebraic equations in Hes- 
senberg form of index vg = 3 (cf. Def. 3.7 and Def. 3.6 for z = [o q A| T and references 


therein). 


In the case of ordinary contact constraints, the constraint forces G74 are ideal-orthogonal 
to the contact constraint manifold C = {q : g(q) = 0} and rank(GM~'G") = d : ns 
applies. 


2 Also frequently referred to as algebro-differential equations, implicit differential equations or singular systems 


((63]). 
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Q GTA 


q(t) 
Tae) Q 


Figure 3.1.: Ideal-orthogonal contact constraint realization. 


The semi discrete equations of motion can then be solved by subsequently integrating 
(3.12) in time by applying suitable finite difference schemes (cf. e.g, [66]). It should be 
kept in mind, that suitable stabilization techniques such as the Gear-Gupta-Leimkuhler 
stabilization (cf. Rem. 3.8) or the Baumgarte stabilization (cf. Rem. 3.9) are mandatory 
for temporal integration with error controlled variable time step size ([10]). 


In Ex. 3.10 and Ex. 3.11, two classical examples of inherent spatially discrete systems that 
align with system (3.12) are presented. 


Definition 3.3 (Constraints). Geometric constraints pertaining the configuration space as 
given by 

g(t,q) =0 (3.15) 
are termed holonomic. Otherwise, the constraints are termed non-holonomic. Examples of 
non-holonomic constraints are, e.g. 


(i) Inequality constraints of the form 


g(t,q) <0. 
(ii) Constraints in differential but non-integrable form 


g(t,q,Drq) =0. (3.16) 


Le. there exists no function G(t,q) such that G(t,q) = g(t, q, Diq). Any integra- 
tion of constraint (3.16) will not lead to a constraint in form of (3.15). An example 
of a constraint in differential representation that is integrable is a non-sliding coin 
on a line (cf. Ex. 3.4). In contrast to that, the same constraint forcing the coin to roll 
without sliding on a plane is non-holonomic (cf. Ex. 3.5). 


Furthermore, a constraint is called rheonomic if the constraint depends explicitly on the time 
variable, i.e. dıg + 0. Otherwise, i.e. 0:g = 0, the constraint is called scleronomic. For more 
details on the classification of constraints see [44, pp. 12-16] and [70, pp. 236-237]. 
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Example 3.4 (Coin on a line: Holonomic constraint). The kinematics of a circular coin 
with radius r rolling on a line without sliding is uniquely defined by the coordinates q = 


[x el" in addition to the scleronomic constraint in differential form as given by 
9(q,Diq) = Dix —r Dig = 0 (3.17) 


Since a function G(q) = x — rọ = const. exists, the constraint (3.17) is holonomic. 


< 


Figure 3.2.: Holonomic constraint: Non-sliding coin on a line. 


x 


Example 3.5 (Coin on a plane: Non-holonomic constraint). Reconsidering Ex. 3.4 for the 
case of a rolling coin on a plane, the kinematics of the coin can be described in terms of its 
location q; = [x y| and its orientation q, = Ç e]. The condition of pure rolling, e.g. 
setting the orientation and the location of the coin into relation, can be ensured by imposing 
the following scleronomic constraints 


(3.18) 


g(q,D:q) = v - r(D;g)n(®) = | vie hee ae 


cos® 


depending on the coordinates q = lq: 4] . If (3.18) would be holonomic, a functional 
relation between q, and q, could be established, hence (3.18) could be uniquely solved by 
integration. This is obiously not the case. Various maneuvers of the coin in a time interval 
T = [0,T] are conceivable such that one coordinate can be chosen independently at time T. 
Consider exemplarily the maneuver of tracing a circle, such that x(0) = x(T), y(0) = y(T) 
and 8(0) = Q(T). Obviously (0) + p(T) depending on the chosen radius of the circle. 


Figure 3.3.: Non-sliding coin on a plane. 
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Definition 3.6 (Differential-algebraic equations in Hessenberg form). Differential-alge- 
braic equations in the form of 


Dizi = Fı(21,..,2,) 
D;Z2 = F2(Z1,..., Zv-1) 
D;Z3 = F3(Za, ES Zy-1) 


D;Zy-1 = F y-1(Zy-2, 15 Zy-1) 
0= F,(zv-ı) 


with a nonsingular 


%,.,F v z Oz, »F y-1 a 0z,F 2 Foca Oz, F 2 
are called DAEs in Hessenberg form of differentiation index va (cf. [63, p.172]). 


Definition 3.7 (Differentiation index). The differentiation index vg of a general differential- 


algebraic equation 
F(z,z,t) =0 (3.19) 


along with a singular ð» F, is defined by the minimum number m € N of derivatives 


D:F(z’,z,t) = (0, F) Dez + (9,)(Dız) + F =0 


D’F(z',z, t) = (dv F) Diz+ ..=0 
(3.20) 


DP F(z’, z, t) = (0, F) Dz +- -0 


that need to be considered, such that Eq. (3.19) can be solved for z'(t) = D;z(t). This 
definition of the differentiation index of general differential-algebraic equations can be traced 
back to C.W. Gear (cf. [40] and [41]). For a comprehensive study of the differentiation index 
see [63, pp.96-113]. 


Remark 3.8 (Gear-Gupta-Leimkuhler stabilization). The original idea of the Gear-Gupta- 
Leimkuhler stabilization is to minimally extend the differential-algebraic system (3.13) by 
temporal derivatives of the constraint equation. These additional constraints are enforced by 
additional Lagrange multipliers. For Hessenberg systems with differentiation index va = 3, 
ie. mechanical systems subjected to holonomic ideal-orthogonal contact constraints, the 
stabilized system in terms of Gear-Gupta-Leimkuhler has differentiation index vg = 2 and 
reads 


Diq=v0-G' pn, 
MD,w=F-G1, 
0=g, 


0=99D:q=Gv. 
This stabilization technique can be traced back to [42]. 
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Remark 3.9 (Baumgarte stabilization). In principle, the second time derivative of the con- 
straint equation 


Dig(g) = 0 (3.21) 


could be used to eliminate the Lagrange multipliers to integrate the differential-algebraic 
system of equation (3.19) with differentiation index vg = 3 numerically. Unfortunately the 
resulting system of ordinary differential equations turns out to be be unstable in the sense 
of Ljapunov and drift-off phenomena may occur. To avoid this, the following stabilization 
method is suggested by J. Baumgarte in [15]. Instead of using only the differentiated con- 
straint equation (3.21), a combination of the original and differentiated constraint equations 
may be taken into account, i.e. 


D’g+2aD,g+ß’g=0, a>0 (3.22) 


in the case of holonomic and 
D;g+yg=0, y>0 (3.23) 


in the case of non-holonomic constraints. In Eq. (3.22) and Eq. (3.23), the non-negative and 
constant parameters a, ß and y have been introduced. 


Example 3.10 (Planar mathematical pendulum). The planar mathematical pendulum is 
an important prototype of a constrained mechanical system. Its virtue lies in encompassing 
the mathematical structure without being distracted by superfluous complexity. A particle 
in a gravitational field g = 1 with mass m = 1 is forced to stay on a circle by a rigid and 
weightless rod with length l = 1 attached to an inertial frame. Its motion is governed either 
by ordinary differential or differential-algebraic equations, depending on whether redundant 
or generalized coordinates are chosen. Both formulations will be addressed subsequently. 


Figure 3.4.: Illustration of a planar mathematical pendulum. 


(i) Redundant coordinates. By choosing cartesian coordinates q = [x yļ”, the La- 
grangian of the planar mathematical pendulum can be established as 


1 1 
L(t, q, D:q) = 3 Diq’ Dıg+y = 5 (Dex) + (Diy)’)+y. 
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Together with the scleronomic holonomic constraint 


ED ae 
9(4) = 5 (4-1) - z +y 1) =0 (3.24) 
and its corresponding Jacobian 
G=|x y]|=g, (3.25) 


the motion of the mathematical pendulum is governed by the Euler-Lagrange equa- 
tions for constrained mechanical systems (cf: Rem. 3.12, Eq. (3.38)) taking the form 
of a system of differential-algebraic equations as given by 


Dix =-Ax, 
D’y=-Ay+1, (3.26) 
g(x,y) =0. 


By introducing 

M=I, F=[0 1] and G=q 
the problem at hand aligns with system (3.12). The system of equations (3.26) is 
sometimes referred to as Lagrange’s equations of the first kind (cf. Rem. 3.12, Ex. 3.13 


and references therein). By introducing the velocity v = D:q a Hessenberg form with 
differentiation index vg = 3 can be identified (cf. Def. 3.6). 


This might be proven by transforming the differential-algebraic system of equations 
into an equivalent system of ordinary differential equations by using derivatives of 
the system and parts of it. In the context of the mathematical pendulum, the differ- 
ential part of system (3.26) might be inserted into the second time derivative of the 
scleronomic holonomic constraint equation (3.24) as given by 


Dég(r(t)) = x D?x + (Dix)? + y Déy + (Diy)? =0. 
This leads to a relation for the Lagrange multiplier given by 
À = (Dix)? + (Diy)? +y. (3.27) 
A further differentiation of Eq. (3.27) with respect to time 
D,A = 2(D,x)(D}x) + 2(Dry) (Diy) + Dry 
= 2A (x Dix + y Diy) + 3 Djy 
=:5 Diy 


gives rise to the following system of ordinary differential equations (cf. [42] and ex- 
ample therein) 


Diqg-v=0, 
D;vo+Ag-F=0, 
D;A-3D,y=0. 


43 


3. Sequential space-time integration 


44 


Thus, the differential-algebraic system of equations (3.26) could be transformed into 
a system of ordinary differential equations by differentiating three times in total. 
Hence, the differential-algebraic system of equations indeed has a differentiation in- 
dex vq = 3. 


In Fig. 3.5, an illustration of an ideal-orthogonal constraint representation is depicted. 
Therein G! A, that enforces the contact constraint, is ideal-orthogonal to the constraint 


manifold C = {r : 9(q) = 3 (q‘q - y) = 0}. 


Figure 3.5.: Ideal-orthogonal constraint representation in context of the mathematical pendulum. 


(ii) Generalized coordinates. Alternatively, the pendulum can be formulated in terms 


of a generalized coordinate p. Therefore, the differential-algebraic system of equa- 
tions (3.26) along with (3.24) can be transformed into one single ordinary differential 
equation. Therefore a suitable coordinate transformation F (Q) : R > R? such that 
q tv F!(o) need to be identified merely. For instance, establishing the transforma- 
tion 


Fy) = | (3.28) 


a relation for the Lagrange multiplier, as given by 


A = cos’ (p) (Drp)? + sin? (p) (Dip)? + cos(p) = (Drg)? + cos(p) 


is obtained. Inserting the transformation (3.28), along with its first and second time 
derivative given by 


D:F(ọ) = a Diy (3.29) 


and 


DEF) = | | Dio + | SNP) | Dio = DFO) - Flo) Die, 


3.1. Direct dynamics problem 


respectively, into equation (3.26) gives rise to 


- sin(p) (D:p)” + cos(y)(D; 9)” = - (Dre)? + cos(p)) sin() . 


Eventually, the system of differential-algebraic equations (3.26) can finally be trans- 
formed into one single ordinary differential equation 


Džo(t) +sin(p) = 0 


governing the motion of the planar pendulum in terms of a single generalized coordi- 
nate. 


Example 3.11 (Overhead crane - direct dynamics problem). Besides the mathematical 
pendulum, analyzed in Ex. 3.10 as a prototype of constrained mechanical systems, the pla- 
nar model of an overhead crane will be considered below (cf. Fig. 3.6 for an illustration 
thereof). This model serves as an important archetype for the inverse dynamics of underac- 
tuated mechanical systems. Therefore, it is first introduced generally in context of the direct 
dynamics problem. 


Figure 3.6.: Planar model of an overhead trolley crane. 


A load with mass mz = 1 is attached to a trolley with mass m; = 1 via an inextensible and 
massless rope. The rope itself is connected to the trolley using a pulley with radius r = 1 
and moment of inertia J = 1, such that the length l : T > R of the rope is adjustable by 
exerting a torque M : T +> R on the pulley. Furthermore, the trolley is capable of being 
maneuvered along a path s € R caused by a force F : T +> R acting horizontally on the 
trolley. The actuation F = [F m|” causes a motion of the load in the x-z-plane. Without 
losing generality, the gravitational constant is assumed to be g = 1. 


Similarly to the mathematical pendulum discussed in Ex. 3.10, the motion of the crane can be 
described either in redundant or generalized coordinates. Both formulations will be addressed 
below. 
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(i) 


(ii) 


f ; ; T 
Redundant coordinates. In case ofchoosing redundant coordinatesq = [s l x z] 


the Lagrangian of the system at hand is obtained as 
1 
L= 5 [ (Drs)? + (Dio)? + (Dix)? + (D:2)°] +z. 


Additionally a scleronomic holonomic constraint given by 


9(q) = HC -s)/'+2-P)=0, (3.30) 


need to be taken into account. Then the Euler-Lagrange equations governing the 
motion of the crane are obtained as 


D?s+A(s—x)-F=0, 

D’I-Al-M=0, 

D’x+Ax-5)=0, (3.31) 
D’z+Az-1=0, 

g9(s,1,x,z) =0 


See Rem. 3.12 and references therein for more details on the Euler-Lagrange equations 
of constrained mechanical systems. The system of differential-algebraic equations 
(3.31) aligns with system (3.12) considering the constraint Jacobian 


G = a49(q) = [-(x-s) -1 (x-s) 2]’ 


as well as F = [F M 0 -1]" and M = I. Consequently, system (3.31) can be 
identified as a differential-algebraic system of equations in Hessenberg form of index 
va = 3. The constraint realization is ideal-orthogonal. 


Generalized coordinates. Instead of describing the motion of the overhead crane in 
redundant coordinates, generalized coordinates might be chosen alternatively. For 
that purpose, a function p(q) can obviously be found, such that x = s + lsin ọ and 
z = lcos»p holds. The motion of the overhead crane sketched in Fig. 3.6 can then 
be reformulated in terms of the generalized coordinates q = [s l gl’. The corre- 
sponding Lagrangian then reads 


L(Diq, q) = - [ (Drs)? + (D;l)? + (D: (s + L sin g))* + (D; (Lcos 9))”| . (3.32) 


In Eq. (3.32) a non-sliding contact of the rope and the pulley, i.e. Dil = rD,o, has 
been assumed. The motion of the crane is then governed by a system of ordinary 
differential equations in the form of (3.11) with 


2 sing Icosp 
M= | sing 2 0 
Icosp 0 1? 
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and 
2(D:1) (Dro) cos p - l? (Dro)? sin F 
F= -l1 (Dro)? - cos -|M 
21(D:1)(D+ọ) - l sin ọ 0 


Remark 3.12 (Constrained mechanical systems). The equations of motion, postulated for 
several (constrained) mechanical systems will be derived in more detail subsequently. The 
variational problem with (finite) constraints based on Hamiltons principle can be stated as 
follows: Find a function q : T > R”, that extremizes the action 


th 
S= / L(t,q,Drq) Dt, (3.33) 
to 


compared to all varied functions 


q =q + ôq =q + an(t), (3.34) 
by additionally complying to holonomi constraints of the form 


g(t,q) =0. (3.35) 


In (3.34), the function dq(t) = an(t) is called the variation of the desired function q(t) 
that extremizes the action (3.33). Thereby, the varied function q(t) is supposed to lie in 
a sufficiently close neighbourhood? of the solution q(t). Thus, the action (3.33) might be 
perceived as a function depending on the coefficient a 


tı 

(a) = / L(t,g(a),D:q(a)) dt (3.36) 
to 

such that (3.36) is extremized for a = 0 compared to all sufficiently small a. Hence, 
Da®(a)laao = 0- 


Concurrently, the constraints (3.35) must not be violated. According to the Lagrange multi- 
plier theorem (cf. e.g. [70, pp. 234-235]) there exist a function A : T > R” such that the 
Euler-Lagrange equations 


-[L]g = DiögL) — gL +A‘ gg =0 (3.37) 


constitute necessary conditions to the function q(t) extremizing the action (3.33). See Ex. 3.13 
for the case ofn = 2 independent functions subjected to holonomic constraints as well as ref- 
erences therein, covering the more general case ofn > 2 functions subjected to non-holonomic 
constraints. Note that the Jacobian G = dgg(t,q) is supposed to be non-singular. Eq. (3.37) 


3 The correlations, that are made subsequently are extendable to the more general case of non-holonomic 
constraints (cf. [48],[27],[46]). 
4 For that purpose the scalar coefficient œ needs to be chosen small enough 
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along with the constraints (3.35) are frequently referred to as ‘Lagrange’s equations of mo- 


tion of first kind’. 


By introducing an augmented Lagrangian in form ofL* = L-A" g, Euler-Lagrange equations 
as given by 
—[L"]q = Di(dgL*) - 9g1* = 0 (3.38) 


can be established. Eq. (3.38) represents necessary conditions to the functions q(t) and X(t) 
extremizing the augmented action 


tı 
Ss“ -f L*(t,q, Diq) dt. 
to 


This is equivalent to (3.37) along with (3.35) (cf. [32], [71] and [68]). 


Example 3.13 (Constrained mechanical systems). Emphasizing the relevance of Rem. 3.12, 
the variational problem with (finite) constraints will be pursued subsequently for the spe- 
cial case of functionals depending on two independent functions subjected to holonomic con- 
straints. The functions x : T œ> Randy: T > R are searched for such that the action 


th 
s= [ L(t,x,y,x’,y’) dt (3.39) 
to 


is extremized, by additionally complying holonomic constraints as given by 


g(t,x,y)=0. (3.40) 


Geometrically, the curves q(t) = [x(t) y(t) |, lying on the constraint manifold C = {q(t) : 
g(q(t)) = 0} and extremizing the action (3.39) are searched for. Throughout this example 
the standard notation f’(t) = D+ f(t) is used to denote the derivative of a function f (t) with 
respect to time t. Assuming non-singular constraint equations, i.e. 0.g + 0 and dyg + 0, 
Eq. (3.40) might be resolved as 

y = k(t, x) (3.41) 


Consequently, the Lagrangian of the problem at hand reads as 
L(t, x,y, x’, y) = L(t, x, k, x’, ðrk + (Ak) x’) . (3.42) 
In Eq. (3.42) the time derivative of (3.41) as given by 
y = D,k(t,x) = &k + (0,k)x’ 


is used. In accordance to classical calculus of variations, a function x(t) that extremizes the 
action 


ti 
S= L(t, x, k, x’, ðrk + (O,k)x’) dt. (3.43) 


to 
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is searched for. Following Rem. 3.12, the first variation of the action (3.43) is obtained as 


ti 
ôS = Da / L(t, Xa, k, x}, rk + (3x k) x!) dt 


to a=0 


th 
= / NOL + (OL) (Ack) + 17 Oe L+ ðwL (3k +x k +y'ðxk) dt 
to 


i (3.44) 
= [el + (WDD) de 
2 t 
+ / (AL + (OL) (Ack) + Oe L (pk + x'a k)) y dt =0. 
Integrating the first integral of Eq. (3.44) by parts gives rise to 
Dy (dL + (3y L) (Ack) — [OL + (OyL) (Ack) + Ay L (xk +x’ Hk) |. (3.45) 
Applying the product rule to the term D: ( (3y L)(Əzk)) and observing that 
Di d,k) = K+ x'a k , 
Eq. (3.45) can be reformulated as 
(D; (aw L) — &L) + (Di (dy L) - dyL) ak = 0. (3.46) 
Since g(t, x, k) = 0 need to be fulfilled at any time, the following applies 
Dug = xg + (dyg)(Axk) = 0. (3.47) 
A comparison of (3.46) and (3.47) inevitably yields 
(D: (ax L) — &L) : (Dy (dy L) - dyL) = 3xg : 3yg . (3.48) 


Since the Jacobian of the constraint is assumed to be non-singular at any time, i.e. 
%g(x,y,t) #0 and dyg(x,y,t) #0VtET, 
there needs to exist a proportionality factor A: T > R such that the necessary conditions 
-[L]x = Dz (AeL) — AL = Adxg and - [L] y =D;(dgL) - ayL =Adyg (3.49) 


to the solutions x(t) and y(t) apply. Eq. (3.49) along with the constraint equation (3.40) are 
also referred to as Lagrange’s equations of motion of first kind. They represent a system of 
differential-algebraic equations, governing the motion of constrained mechanical systems. 


These equations of motion are equivalent to the Euler-Lagrange equations 
- [L"], = Dae") - a,L* = 0, 
—[L"], = Day L") — ayL"=0 and 
- [L*], =—-a,L* = -g(t, x, y) = 0, 
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representing necessary conditions to the functions x(t) and y(t) as well as to the Lagrange 
multiplier A(t) extremizing the action 


ti 
S* =f L* (t,x, y, x,y) dt (3.50) 


to 


for the augmented Lagrangian L* = L-Ag For more details on variational problems subjected 
to (finite) constraints, see [32, pp.143-220], [71] and [68]. 


3.2. Inverse dynamics problem 


While the classical direct dynamics problem aims to find the overall motion q(t) of a 
(constrained) mechanical system excited by a given actuation f(t), the inverse dynamics 
problem seeks for an actuation f(t) such that the motion of the mechanical system can 
partly be specified. 


In this context, inverse dynamics problems might be interpreted as mechanical systems 
subjected to rheonomic holonomic constraints. Neglecting any further geometric con- 
straints, the motion of such systems is governed by differential-algebraic equations in 


form of 
M D;q + F(D:q.q,t) + By f(t) =0, 


h(t,q) = Hq - y(t) =0. 


In (3.51)2, the matrix H is of boolean type. Its purpose lies in extracting the prescribed 
degrees of freedom from the nodal configuration vector q : T +> R*. The actuating com- 
ponents f : T > R” then take the role of Lagrange multipliers enforcing the imposed 
servo-constraint? conditions (3.51)2. Note that m < k is assumed, since we are focusing 
on underactuated systems. 


(3.51) 


In contrast to (ideal)-orthogonal contact-constraints that have been discussed in Sec. 3.1, 
servo-constraints in general do not have collocation property. Geometrically, this means 
that the Lagrange multipliers enforcing the constraints are not orthogonal to the con- 
straint manifold anymore (cf. [22, 87]). 


We will demonstrate, that the spatial disjunction of the constraint and the constraint 
forces, i.e. the non-standard construction of the constraint realization, causes differential- 
algebraic equations, that are characterized by either high differentiation index’ or the 
appearance of (unstable) internal dynamics. 


> Servo-constraints are also frequently referred to as control- or program-constraints, (cf. [79], [64] and [22], 
respectively.) 

6 E.g. concerning the imposition of essential boundary conditions in the finite element methods, the constraints 
are commonly assumed to have collocation property [3]. 
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Although, servo-constraints could be successfully applied to low-dimensional underac- 
tuated mechanical systems such as cranes [4, 26, 75, 98] and manipulators with passive 
joints [14, 18, 23], the treatment of higher dimensional systems is still an open issue. 


In particular, we face the problem that both, (unstable) internal dynamics as well as high 
differentiation index’ along with high demands on the smoothness of the prescribed tra- 
jectory are affected adversely by the spatial discretization. This and its impact on the 
solvability of the problem at hand will be elaborated subsequently. 


Following [24], the geometrical properties ofthe constraint realization might be specified 
in terms of 
n= rank(HM"' By) (3.52) 


such that three distinct cases of the orientation of the actuation B; f(t) on the constraint 


manifold C = {q : g(q, t) = 0} can be identified. 


For n = m, the constraint realization is referred to as (non-ideal) orthogonal. In conse- 
quence, this leads to differentially non-flat systems where (unstable) internal dynamics 
may arise. View Fig. 3.8 illustrating non-ideal orthogonal constraint realizations in con- 
trast to ideal orthogonal constraint realizations, depicted in Fig. 3.7. (cf. [56] or [26]). The 
notion of differential flatness (cf. Def. 3.23 and references therein) goes back to [38] while 
the concept thereof can be traced back to [49]. 


B; f 
C (dgh)A 


q(t) 
Tq(t)C 


Figure 3.7.: Ideal orthogonal constraint representation. 


Forn = 0 or 0 < n < m, the constraint realization is called tangential or mixed orthogonal- 
tangential, respectively. In the case of tangential constraint representations, the actuation 
Bf is tangential to the control constraint manifold C = {q : g(q) = 0}, i.e. none of the 
constraint components can be actuated directly. In contrast to that, n constraint compo- 
nents can be actuated directly in the case of a mixed orthogonal-tangential representation. 
In both cases, the system at hand is possibly differentially flat without internal dynamics 
or non-flat with internal dynamics See Fig. 3.9 for a geometrical interpretation of the tan- 
gential constraint representation. In Tab. 3.1, the geometry of the constraint realization 
depending on Eq. (3.52) is summarized. 
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(dgh)A 


Cl 


Figure 3.8.: Non-ideal orthogonal constraint representation. 


In the case of differentially non-flat systems, (unstable) internal dynamics may occur, 
hindering numerical integration of the problem at hand. Therefore, it is inevitable to 
identify any (unstable) internal dynamics and carry out relevant analysis of the resulting 
(non-) minimum phase systems (cf. [17] and [85]). On the other hand, flat systems lead 
to differential-algebraic systems that are characterized by a high differentiation index vq 
(cf. Def. 3.7 and references therein). The concept of characterizing differential-algebraic 
systems of equations by its differentiation index vg is highly related to the concept of 
relative degree r (cf. Def. 3.14 and references therein). 


B; f 
C (dgh)A 


q(t) 


Tg(t)C 


Figure 3.9.: Tangential constraint representation. 


Since a numerically stable solution to the resulting differential-algebraic system of equa- 
tions is depending significantly on the differentiation index, the differentiation index 
need to be reduced in order to achieve stable numerical solutions (cf. [4] and references 
therein). 


(i) n = m: (non-ideal) orthogonal constraint realization, ie. all m con- 
straint components can be actuated 


(ii) 0 < n < m: mixed orthogonal-tangential constraint realization, i.e. 
only n constraint components can be actuated directly 


(iii) n = 0: tangential constraint realization, i.e. none of the constraint 
components can be actuated directly 


Table 3.1.: Geometry of the constraint realization depending on n = rank( HM Br) 
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Subsequent examples will demonstrate that both internal dynamics and high differentia- 
tion index are a direct consequence of the discretization process. This obviously restricts 
the applicability of the classical semi-discretization approaches. 


Definition 3.14 (Relative degree). Consider a single input - single output system as given 
by 
D;z(t) = F + Bu(t) , 


(3.53) 
y(t) = H(z) . 
System (3.53) is said to have relative degree r at Zo if 
LgLFH(z) =0, LglZH(z)=0, ++: Lgl;H(z)=0Vi<r-1 (3.54) 


along with a non-singular Lele H(z) # 0 applies. In Eq. (3.54), 


LFH(z) = SH) 


a=0 


denotes the Lie-derivative of H along F (cf. e.g. [56, Sec.1.2.]). The relative degree r of (3.53) 
can be interpreted as the number of derivatives of the output y(t) that has to be taken such 
that the input F(t) appears explicitly (cf. [56, Sec.4.1. p.139]). With (3.54), a coordinate 
transformation can be established such that system (3.53) can be transformed into a Byrnes- 
Isidori input-output normal form (cf. [86, Sec.3.1.]). 


Example 3.15 (Spring-mass system mounted on a moveable carriage). The following 
example aims to underline the importance of the geometry of the constraint realizations. 
Therefore, we aim to analyze the following control problem: On a horizontally moveable 
carriage with mass m; = 1, a harmonic oscillator with mass ma = 1 and stiffness k = 1, 
rotated relatively to the carriage by an angle ọ, is attached. Furthermore, the external load 
F : TR is assumed to act horizontally on the carriage. In Fig. 3.10, an illustration of the 
model is depicted. 


Figure 3.10.: Illustration of a spring-mass system mounted on a moveable carriage. 
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Consulting Rem. 3.12 the Lagrangian 


N 


er 

U U Ss 

L=++24- > (3.55) 
2 2 2 

gives rise to the Euler-Lagrange equations given by 


2 D’xı + cos(@) Dis -F=0 and cos(p) Dex, + Dis +s=0. (3.56) 


In Eq. (3.55), the velocities of the carriage and the mounted weight have been introduced as 


Enz ea TE ee k +cos(p) Dis (3.57) 


0 sin(p) Des 


respectively. Regarding the inverse dynamics problem, the actuating force F(t) is searched 
for, such that the horizontal motion x2 of mass mz follows a prescribed trajectory as given 


by 
h(x, s,t) = x2 — y(t) = x1 +scos(p) - y(t) =0. 


Analyzing the geometry of the constraint realization, three distinct cases can be observed: 
(i) @ = 0: Ideal-orthogonal constraint realization - stable internal dynamics occur 
(ii) 0 < g < 3: Non-ideal orthogonal constraint realization 
(iii) ọ = 5: Tangential constraint realization - no internal dynamics occur 
Further analysis regarding this control problem can be found in [23] and [87]. 


Example 3.16 (Linear elastic bar - inverse dynamics problem). Revisiting the linear- 
elastic bar, introduced in Ex. 2.5 and Ex. 3.1 for the pure Neumann problem, the semi-discrete 
equations of motion (3.4) in conjunction with servo-constraints boil down to a system of 
differential-algebraic equations given by 


Mij Diu; + Kjjuj = —oinf A 


Up+ ZY. 


(3.58) 


Subsequently, a consistent approximation is compared to a lumped approximation of the 
inertia terms with regard to its impact on the properties of the constraint realization. 


(i) Lumped mass matrix. Consider a mass matrix with entries Mi; = Miôij. Regarding 
the last two nodes p and p + 1 associated with the p-th element, the differential- 
algebraic system of equations (3.58) yields 

Mp Diup + Kp pup + Kppriy = 0 Er 

Mp+ı Diy + Kpsi,pulp + Kpsipuiy = 0. 


Here, the summation convention does not apply. Since y(t) is prescribed, it can be 
deduced from the last equation that the displacement of node p is given by 


Up = -Korp (Mp+ı Diy + Kp+1,p+1Y) = Up (Diyy) . (3.60) 
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Consequently, 

Drup = ap (Diy, Diy) (3.61) 
is obtained where the function ap follows directly from differentiating (3.60) twice 
with respect to time. Repeating this procedure for the elements p—1,...,p—n reveals 
that 


Up-n = Up-n (Dey; D?"y, o’ Y) Deus = Ap-n (Des iD rae ree D?y) 
Setting n = p — 1 leads to 
u = u (DnD Y, Late r) : Diu; =a, PYDI D Dir) (3.62) 
so that for i = 1, (3.58); yields 


f ID YD Dr r) ; (3.63) 


Differentiating Eq. (3.63) once again with respect to time, reveals that the time deriva- 
tive of the algebraic variable f pertaining to the differential-algebraic system of equa- 
tions (3.58) can be obtained after 2p + 3 time derivatives of the algebraic constraint 
(3.58)2. This implies that the system (3.58) has differentiation index vg = 2p + 3 
(cf. [13, 63]). Accordingly, the index increases with the number p of finite elements 
used for the space discretization of the bar. Since all of the unknowns of system (3.58) 
can be expressed in terms of y and derivatives thereof, the present problem exhibits 
the property of differential flatness ([67, 90]). In particular, the nodal displacement 
Up+i plays the role of a flat output. Analogous results have been obtained for the ar- 
ticulated mass point systems dealt with in [22, 26, 74]. Note that the simple result for 
the flat output in Ex. 3.16 hinges on the mass matrix being diagonal. The situation 
gets more intricate in case of a consistent mass matrix or higher-order finite elements. 


Special case n = 1. Using only one finite element, the problem boils down to the 
system of ordinary differential equations 


Dĉu +u -u =-f, 


Dîu - u +w =0. 
subjected to a rheonomic holonomic constraint given by 
gu) =u — y(t) =0. 


Essentially, this differential-algebraic system of equations governs the inverse dynam- 
ics of a vibrating (double-) oscillator. By introducing 


Ie ofl wef 


55 


3. Sequential space-time integration 


as well as the actuating force f acting as Lagrange multiplier, the resulting system 
aligns with system (3.51). Due to the non-collocation of the constraint force and the 
imposed constraint, the constraint realization is referred to as tangential. In the case of 
this two-dimensional system, the constraint realization can be visualized in Fig. 3.11. 


ug A 


HA 


B; f 
YOH 727777777” 


u 


Figure 3.11.: Tangential realization of the servo-constraint. 


(ii) Consistent mass matrix. For the more general case of a consistent mass matrix (cf. 
Ex. 3.1) and confining to piecewise linear Lagrangian shape functions L;(s) it can be 
observed thatduetoH=|0 0 0 >> 0 1| and Br = [1 0 0 > 0 0| 
the norm of P = HM"'B, depends on the number of elements p. Apparently, P 
becomes singular as p tends towards oo. Hence, by increasing the number of elements, 
the constraint realization tends to become tangential. (cf. Fig. 3.12) 


1.5 


Pl 


1 4 7 10 
P 


Figure 3.12.: Norm of P = HMB; depending on the number of elements p. 


Example 3.17 (Overhead crane - inverse dynamics problem). Revisiting the model of the 
planar overhead crane introduced in Ex. 3.11, the corresponding inverse dynamics problem 
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aims to find an actuation f = [ F(t) Ma] such that the attached load follows a pre- 
scribed trajectory 


x(t) = yx(t) , 
z(t) = yz(t) . 


Similar to the formulation of the direct problem, the inverse dynamics can be established 
either in redundant or in generalized coordinates. 


(i) Redundant coordinates. Using redundant coordinates q = [s l x 2)", the control 


(ii) 


problem can be formulated in terms of a differential-algebraic system of equations as 
given by 
Diq+ F—G'A-Brf =0, 


g(q) = e-t -P)=0, (3.64) 
h(q) =0 


where F = [0 00 |" denotes the external applied excitation, excluding the 
actuation f. Additionally to the geometric, i.e. ideal-orthogonal constraint g(q), that 
is enforced by the constraint forces GA = (dgg)A (cf. Ex. 3.11), the system at hand is 
subjected to servo-constraints 


or] so 


that are enforced by the external actuation Bff. Due to the non-collocation of 


1000 
Bl to ol 


and the Jacobian of the servo-constraint 


(3.66) 


aha) == |) 0 1 | 


0 0 0 1 


it applies, that n = rank(HM"! BT.) = 0. Hence the realization of the servo-constraint 
can be identified to be tangential (cf. Fig. 3.9). As a consequence, the differentiation 
index of the underlying differential-algebraic system of equations can be determined 
to be vg = 5. 


N ; ; ; : T 
Generalized coordinates. In the case of using generalized coordinates q = [s l o] ; 


the inverse dynamics problem is given by the differential-algebraic system of equa- 
tions 5 P 
MD;q=F-B;f, 


(3.67) 
h(q) =0 
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for 


_|x-Yyx} _ |stlsin(g) — yx} _ 
h(q) = i 7 >| =I Bes) |” 0 (3.68) 
and 
1 0 0 
Br = f i T (3.69) 
Again, the non-collocation of Bf and 
7. |1 sin(g) lcos(ọ) 
= b cos(p) -Isin(p) 


leads to rank(HM~'B;) = 1 and consequently to a mixed orthogonal-tangential 
constraint realization. 
Remark 3.18 (Index reduction by minimal extension). By augmenting the differential- 
algebraic system at hand by feasible index 1 systems the differentiation index can be reduced. 
Since the resulting system is singular so-called dummy variables, i.e. additional unknown 
functions, need to be introduced. This concept of index reduction by minimal extension goes 
back to [73] and is motivated by the GGL stabilization [42]. 
Example 3.19. Consider a Hessenberg system with index vg = 3, taken from [73], as given 
by 
D;u(t) = v(t) , 
Dolt) = f(t), (3.70) 
g(u, t) = u(t) — y(t) =0. 
This system might be augmented by successive derivatives of (3.70), and (3.70)3, i.e. by a 
differential chain in index 1 form given by 


D;u(t) = Dey(t) , 
Dju(t) = Dey(t) , (3.71) 
D;,o(t) = Déu(t) . 
System. (3.70) along with system (3.71) form an overdetermined system for the unknown 
functions u(t), v(t) and f(t). Therefore, dummy variables & = D;u(t), & = D?u(t) and 


& = D,v(t) need to be introduced in order to be able to solve the resulting system of purely 
algebraic equations (vq = 1) uniquely. 


Example 3.20 (Bead on parabola). The motion of a bead sliding on a parabola is governed 
by a system of differential-algebraic equations given by 


Du — u = 0, 
Druz -02 =0, 
Divi — 2Au =0, (3.72) 


D;v2 -14+A=0, 


glu, U2) = uz -u$ =0. 
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System (3.72) is characterized by a differentiation index vg = 3. Augmenting system (3.72) 
by the first and second time derivative of the constraint equation (3.72); with respect to the 
time variable, i.e. 
D;g(u1,u2) = v2 — 2vıu = 0 
and 
D’g(u,, u2) = Dvs - 2 D,vıuı — 2D,07 =-(1+/A)- 4u? — 2v; = 0, 
respectively, it can be transformed into a fully algebraic system of equations, i.e. vg = 1, by 


introducing the new variables & = D;uz and & = D;v; (cf. [63, Sec. 6.4]). 
Example 3.21 (Mass-spring system - index reduction by minimal extension). Revisiting 
Ex. 3.16 and considering the special case of using one finite element for the spatial approxi- 
mation of the solution, yields 
Du!) + u(t) — uz(t)- f(t) =0, 
D?up(t) — u(t) + u2(t) =0, (3.73) 
ga) =u — y(t) = 0. 
The system of differential-algebraic equations (3.73) features a differentiation index vq = 5. 
Augmenting system (3.73) by D?u2(t) = Diy(t) = &(t) yields 
Deus (t) + u(t) — un(t) - f(t) =0, 
č (t) — u(t) + ua(t) =0, 
glu) = uz — y(t) =0, 
Alt) = Dy). 
A Hessenberg form with differentiation index va = 3 can be identified by introducing the 
velocity vı = D:u, (cf. Def. 3.6). Augmenting system (3.74) by 
Diu (t) = Diy) + y(t) = & 0) (3.75) 
a fully algebraic system of equations featuring va = 1 as given by 


&2(t) +u (t) — u2(t) — f(t) =0, 
& (t) —ui(t) + u2(t) =0, 
gluz) =u —y(t) =0, (3.76) 
E(t) =Diy(t), 
&(t) = Dfy(t) + y(t) 


can be established. By using the constraint equation (3.73)3, the prescribed function uz serves 
as a parametrization, to transform the differential part of system (3.73) into a purely alge- 
braic system of equations and finally solve the problem at hand without actually integrating. 
In this context, the parameterization uz serves as a differentially flat output algebraizing the 
differential-algebraic system of equations (cf. Def. 3.23 as well as accompanying Ex. 3.24, 
Ex. 3.25 and Ex. 3.26 for more details on differential flatness). This indicates a conceivable 
minimal extension of the equations, aiming to reduce the differentiation index. 


(3.74) 
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Example 3.22 (Overhead crane - index reduction by minimal extension). Revisiting the 
overhead crane model introduced in Ex. 3.11 and Ex. 3.17 in the context of the direct 
and inverse dynamics problem, respectively. In case of choosing redundant coordinates 


q= [s l x z| a the actuation f(t) satisfying the system of differential-algebraic equa- 
tions given by 
Dig - ögg(q)A- Bff = 0 

gg) =0 (3.77) 

h(q) =0 
is searched for. Augmenting system (3.77) by the second time derivative of the servo-constraint 
given by 

D;h(q) = 0 


and introducing the new variables & = D?x and & = D?z, yields a system with differ- 
entiation index vq = 3. An algebraic relation for the Lagrange multiplier A is then given 


by 
1 
A=—(1-&). (3.78) 
Ya 
Definition 3.23 (Differential flatness, cf. [49]). A differential equation of the form 
F(t, u(t), Dyu,--- , DP u, f(t), Dif,- , DP f) =0 (3.79) 


is called differentially flat, if there exists an arbitrary function € € R + w(é) € R with 
derivatives Drw(£), +++ ‚Diw(8) such that a transformation of the form 


u = a(é,w, Dew, Dew, <) 
f = b(& w, Dew, Dew, +) (3.80) 
t = c(& w,Dew, Dew, +) 


satisfies the differential equation (3.79) by taking into account 


D = _ Dza 2 5 
u(t) = Dza Deč = Dic’ D;u(t) = ——(D; za Dec - Dza Dc), 


T: c)? 


Deb 
D:f(t)=DzbD.é=—-, Dif(t Dzb Dec - Deb Dic), 
if (t) =DebD.E De fO = D: Day ec — Deb Déc) 
along with 
Dza(é, w, Dew, EN +) = Ira + ðwaDgw + (Dew) a Dew + 
Dzb(é, w, Dew, Dew ++) = ðzb + dy b Dew + Anew) b De wt: 


Alternative definitions of differential flatness can be found e.g. in [67] and [38]. See also [90] 
for more background on linear systems. 
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Example 3.24 (Differential flatness). Consider the first-order ordinary differential equa- 
tion 
F(t, Diu, f) = Du(t) - f(t) =0. (3.81) 


The transformation 


u = a(é,w, Dew, Dew, --)=w(E), 
f = b(č, w, Dew, Dew +) = D:w(d),, 
t = c(& w, Dew, Dew )=E 
with 
D;u(t) = Dew(¥) 
satisfies the differential equation (3.81). 


Example 3.25 (Differential flatness). Consider the second-order ordinary differential equa- 
tion 
F(t, D?u, Dyu, f) = D?u(t) + u(t) - f(t) =0. (3.82) 


The transformation 
u = a(ğ, w) = w(d) 
f =b(¥, w, Dew, Dw) = Diw(d) + w(€é) 
t=c(d)=8 
with 
D;u(t) =Dew(£) and Diu(t) = Dyw(d) 
satisfies the differential equation (3.82). 


Example 3.26 (Differential flatness). Consider the system of second-order ordinary differ- 
ential equations 


Duy + ui — U2 =f, 
se (3.83) 
Di uz + uz - u =0 


governing the motion of a mass-spring system analyzed in Ex. 3.16 and Ex. 3.21. The trans- 


formation 
uy = a(ğ, w, Dew) = Dew(€) + w(Z), 
uz = b(E,w(Z)) = wi), (3.84) 
t=c(%) =% 


allows solving the differential equation without actually integrating it. Inserting the trans- 
formation (3.84) satisfying (3.83), into (3.83), yields a relation for the function f : T >» R 
given by 

f = d(& DIw(d, DEw(H)) = Diw(d) +2D!w(d). (3.85) 
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Remark 3.27 (A spectral Galerkin approach). Alternatively, for interpolating the test and 
trial functions (3.3) in the weak formulation (3.2), global trigonometric basis functions ofthe 
form 


2 : f 
L(s)=1, Ls(s) = = - 27 +1 and L;(s) = sin( 5) — G -= sL) 


forj=1,:-- ,p can be used.’ See Fig. 3.13 for an illustration of these basis functions. Subse- 
quently it will be demonstated that indeed high differentiation index might be circumvented 
but eventually (unstable) internal dynamics occur. 


Lj 


Figure 3.13.: 
Global trigonomteric basis functions for p = 5 


The resulting differential-algebraic system of equations can then be obtained as 


M;i; Did; + Kijdj - B; f = 0, 


(3.86) 
h= y(t) - Hid; =0. 
In Eq. (3.86), the two Boolean matrices Bf and H are defined as 
Bp=[I I 0 --- 0| and H=[I 0 0 > o], (3.87) 


respectively. Furthermore Kı; = Kj, = 0 applies. If n = rank(HM~'B;) = m, where m 
is the number of the actuating components, the constraint realization is called (non-ideal) 
orthogonal (cf. Fig. 3.14 and Fig. 3.15 illustrating the difference of ideal- and non-ideal- 
orthogonal constraint realization). Orthogonal constraint realizations generally lead to dif- 
ferentially non-flat systems where (unstable) internal dynamics may arise? (cf. [56] or [26]). 


7 The idea of using this trigonometric ansatz functions has been suggested by R. Altmann after the talk [94]. 

8 Forn = 0 or 0 < n < m, the constraint realization is called tangential or mixed orthogonal-tangential 
respectively. In both cases, the system at hand can be differentially flat (without internal dynamics) or non- 
flat (with internal dynamics). 
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dz 
HTA 


T 
Br 


V(t) AAA APPT TTT TTT 


dı 


Figure 3.14.: Non-ideal- orthogonal constraint realization. 


Due to the positive definiteness of the mass matrix M for (3.86) it applies thatn = m. Thus the 
control is non-ideal orthogonal to the constraint mainfold (tangential components exists, cf. 
(3.88)). Hence, the differentiation index vg = 3. Therefore the output cannot be differentially 
flat and (unstable) internal dynamics may occur. 


dy À 


y(t 


dı 


Figure 3.15.: Ideal- orthogonal constraint realization. 


The internal dynamics of the system (3.86) may be identified by taking the second time 
derivative of the algebraic constraint (3.86)2. Inserting it into the differential part (3.86), the 
actuation f(t) can be established as 


f=Mu D’y +My Did; +My Did; Perea (3.88) 


A dependency of f(t) on the unknown functions d;(t) fori € {2,3,...} can be observed. By 
inserting the relation (3.88) into the remaining differential part of (3.86), the solution d;(t) 
is governed by a system of ordinary differential equations given by 


M2-Mı M23-My3_ ---| [Diaz Kz Ka ---| | do (Moi - Mi)y(t) 
M3—-Mı2 M3-Mıs -| |D?d; 4 |K32 K33 +| |d3| — | (M31 - Mı1)y(t) 


The system (3.89) constitutes the internal dynamics of the control problem at hand (cf. [26]). 
Further investigations, regarding the stability of the internal dynamics (3.89) are mandatory. 
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The existence of internal dynamics indicates that the present output (i.e. the displacement 
of the elastic rod ats = L) cannot be a flat output of the discrete system at hand. This 
observation leads to the question whether there does exist a flat output. That is, given the 
discrete system with actuating force at the left end, does there exist an output which turns 
the system differentially flat. This might be figured out by introducing v(t) = D;d(t). Then 
the equation of motion subjected to the servo-constraint can be transformed into a first-order 
system in state-space representation 


Dıx = Ax + bf (3.90) 


with 


fh bale = ola 


Following [90], a state-space coordinate transformation of the form z = Tx can be carried 
out by using the inverse of Kalmans controllability matrix 


T=|b Ab A&b |". (3.91) 
This yields a system of ordinary differential equations in new coordinates as given by 
Diz = Az+Tbf . (3.92) 


For the transformed system, it is possible to find the differentially flat output corresponding 
to the actuation f (cf. Ex. 3.16 for p = 2 and Ex. 3.29 for p = 3). 


The proposed global trigonometric discretization of the weak formulation leads to a differ- 
ential-algebraic system characterized by a differentiation index vg = 3. However, the ac- 
tuating force f(t) cannot be identified as a differentially flat output corresponding to the 
partly specified motion y(t). Consequently, internal dynamics that are eventually unstable 
occur. Further investigations should be pursued to elucidate the implications the internal 
dynamics may cause. Incidentally, a discretization of the weak formulation at hand by us- 
ing Lagrangian shape functions assuming a consistent mass matrix also leads, at least for 
small p, to differential-algebraic equations with differentiation index vg = 3 and internal 
dynamics (cf. Ex. 3.28). 


Example 3.28 (Linear elastic rod - trigonometric discretization for p = 2). Consider a 
linear elastic rod with length L = 1, axial stiffness EA = 1 and mass density pA = 1. By 
using p = 2 ansatz functions, the coefficients of the semi-discrete equations of motion can be 


identified as 


|1 1 _ | 0 0 _|1 _ |} 
M= |e AE KE sn, vk sr =|; and ft ||: 


Due to 
n= rank(HM_~'B;) =1 
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the actuator is (non-ideal) orthogonal to the constraint manifold (cf. Fig. 3.11). Differentiat- 
ing the algebraic constraint twice and inserting it into the differential part of Eq. (3.86) the 
actuating component is obtained as 


f(t) = = Did, (0) +Déy(t) . (3.93) 


Taking the time derivative of (3.93), an ordinary differential equation for the actuating force 
f(t) depending on d2(t) can be established as 


D: f(t) = = Did, (0) +D3y(t) . (3.94) 


By inserting (3.93) into the remaining differential part of (3.86), it becomes apparent, that 
the concerned function d,(t) is governed by a ordinary differential equation given by 


1 
r Did; — dz =—-D?y(t) . (3.95) 


Eq. (3.95) can be identified as the internal dynamics of the problem at hand. In this case the 
internal dynamics are obviously unstable. Subsequently we will try to find the differentially 
flat output corresponding to the actuation f (t). This may lead to a better understanding of 
the problem at hand and its controllabillity conditions. Following Eq. (3.90) - Eq. (3.92), a 
transformation of the state-space variables can be obtained by establishing Kalmans control- 
lability matrix as 


01 0 -55 000 0 1 
01 0 35 100 0 0 

-1_ = — 

T =j o -s5 op A’lbıos and Tb=|, 
1 0 15 0 001 0 0 


The actuating force f(t) can then be determined as 
ft) = Diz, - 15 D}z.. 
The transformation z = Tx yields z4(t) = 7 (dp —d,), hence 
f =f (Did), D?d,, Did}, D2d,) . 


The function z4(t) = 4 (d: — dı) represents a differentially flat output corresponding to the 
actuation f(t). Consequently, in case of not choosing z4(t) as output, (unstable) internal 
dynamics may occur (cf. Eq. (3.89)). 


Example 3.29 (Linear elastic rod, p = 3). By using p = 3 ansatz functions, the coefficients 
of the semi-discrete equations of motion are obtained as 


1 1/3 2/r — 7 f6 0 0 0 1 
M= Ys m +! —*/ol, K= 4h "fa An and B’=|ı 
Sym. 1/2 — 8/q2 + 7° fo Sym. 5x” je — 8 0 
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Again, the state-space variables can be transformed using Kalmans controllability matrix 


and Tb= 


Oro 0oo0 
coo cor Fe 
OSS: OOS 
OO O OG et 8 

au oo G 
oo, © 
coo oro co 
oOOHr 000 
Oor. ooo oo 
=. oooo oo 
cocoon oo 
O oOo ooom 


f 


where a, b,c,d,e, f and h are constant values. A dependency of the actuating force f(t) on 
Z6 = ze (1, do,d3) and its time derivatives can be observed, i.e. 


f = Diz + D’z, = D’z, 


Transforming the coordinates back into the original state-space variables indicates that a 
linear combination of dı, dz and d3 represents a differentially flat output. 


f = f(D/d,,D!d;, D!d;, D}d,, D}d;, D}d;, D?d, , D?dz, D?d;) 


A full paramterization of the actuating force f (t) by only using dı (t) and its time derivatives 
is likewise not possible. 
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integration 


Due to the highly restrictive applicability of solving the inverse dynamics problem se- 
quentially in space and time, this chapter aims to analyze the initial boundary value 
problem in more detail. In particular, by exposing the underlying hyperbolic structure 
ofthe governing partial differential equations it is anticipated to gain more insights into 
the problem at hand. Enlighting resulting mechanisms of the flow of informations within 
continuous structures will pave the way to a numerically stable integration of the inverse 
dynamics problem at hand. In Sec. 4.1 therfore, the initial boundary value problem that 
has been postulated in Sec. 1.1, is analyzed in more detail. The new insights will motivate 
to develop in Sec. 4.1 and Sec. 4.2 novel numerical methods that are based on a simultane- 
ous space-time discretization of the governing equations. In Sec. 4.3 numerical examples, 
including mechanical systems, that already have been introduced in Sec. 2.1, Sec. 2.2 and 
Sec. 2.3, are presented underpinning the relevance of the proposed methods. 


4.1. Wave propagation 


Spatially 1-dimensional systems. Considering 1-dimensional systems, the wave equation 
for the control problem (1.1), that has been introduced in Sec. 1.1, can be transformed 
into a system of first-order partial differential equations 


Adv — ð (Bp) =C, 


(4.1) 
Bo; p = Bð, v =0 


by defining the functions v(s,t) = ox(s,t) and p(s,t) = 9x(s,t). Making use of the 
product rule, ie. Bo; p = 0;(Bp) — 9,Bp, Eq. (4.1), is obtained as 


(Bp) = Bo,v = 9Bp . (4.2) 
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By taking advantage of the equality of mixed partials!, hence 0;p = 0,v as well as using 
the chain rule for 


0; B(p(s,t)) = (3p @ B)- a p= grad, (B) "IP, 
the following relation for (4.1), can be established 
9, (Bp) — Bd,v = (dp @ B)posv. (4.3) 


Then Eq. (4.1)2 and Eq. (4.3) are forming a system of first-order partial differential equa- 


tions given by 
A 0||ov 0 I||vo| _|C 
o z||Bp|, |H 0||Bp|, [0] ve) 


In Eq. (4.4), the function H : Q > R@@ has been defined as 
H(p)=B+(, ®B)p. (4.5) 


By introducing the column vectors z = [v Bp]! : Q =œ R”! and F : Q œ> R? as well as 
the square matrices D : Q œ R?%2d and E : Q m R?%2d, Eq, (4.4) might be written more 
compactly as 


Do;z+ Ed,z =F. (4.6) 


Assume there exists a line s = k(t) along that a solution z = z(k(t),t) = Zo(t) is given. 

Then this line is called a characteristic line if the partial derivatives of the solution cannot 

be uniquely determined through given informations along this line, ie. Eq. (4.6) along 
with 

k(t) = 2 4.7 

IZ FAZ T (t) = qt) , (4.7) 


hence 
E D K(t OsZ I D Zo t 4.8 


cannot be solved uniquely for the partial derivatives 0,z and ð;z. If the determinant of 
the coefficient matrix 


d 
Q=E- Da k(t) (4.9) 


as well as the determinants of the matrices, defined by exchanging columns by the right- 
hand side of Eq. (4.8), vanish then according to Cramer’s rule, Eq. (4.6) can be transformed 
into a system of ordinary differential equations along characteristic lines, viz. 


det(Q) = 0 and det(Q) =0, (4.10) 
j 


1 also frequently referred to as Schwarz’ theorem (see e.g. [8, Chap. 24.1 pp. 879-881]). 
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where det;(Q) denotes the determinant of Q where the the j-th column of © is replaced 
by the right hand-side of Eq. (4.8), e.g. the components of the j-th column of Q are given 


by Q; = (F = D#z0(1)) . 


Eq. (4.10); and Eq. (4.10)2 are referred to as directionality and compatibility condition, 
respectively. The procedure outlined above, is also referred to as Method of Characteristics 
([1, 31, 33, 57, 84]). The resulting system of ordinary differential equations can then be 
solved either (if possible) analytically or graphically-numerically (cf. Rem. 4.8). 


For a more comprehensive geometrical interpretation of the underlying structure, see 
Rem. 4.12 and Rem. 4.13 as well as the accompanying examples regarding both (quasi-) 
linear (Ex. 4.14) and general non-linear (Ex. 4.15 and Ex. 4.16) partial differential equa- 
tions. 


Spatially a-dimensional systems. In principle the same holds true for spatially higher 
dimensional systems, i.e. spatial variables s € S c R” for æ > 1. Introducing P = B - 0.x, 
the second term of the partial differential equation (2.11), introduced in Chap. 2 is the 
divergence of P(s, t) as given by 


divs P= div,(B . Osx) = Os, (Bij9s,X;) = Os, (Bj js, x;) rief Os, (Bij 9s, xj) $ (4.11) 
Then Eq. (1.1) takes the form 
Ad; x — (ds, (Bas x) + +++ + ds, (Bas,x)) =C. (4.12) 


By introducing v(s,t) = ôx and p(s,t) = 0.x and making use of the equality of mixed 
partials, i.e. 0:p = 9sv the partial differential equation (4.12) can be transformed into a 
system of first-order partial differential equations given by 


Dð;z - (ds, (E12) +++++9s,(Enz)) =F. (4.13) 


In analogy to the spatially one-dimensional case, characteristics can be found, by search- 
ing manifolds in the space-time domain on which initial data leads to an indifferent solu- 


tion z = [o pl T of the initial boundary value problem. Consequently the inner deriva- 


tives 
Os,Z = Ds2- (AZ) (Ds, t), 


Os, Z = D,„z T (9:2) (Dsn t) 


give rise to the following system of equations 


pele | | en pace 
a a ds; dsn 
By evaluating the determinants of 
Sı dsn 
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the system of partial differential equations (4.13) can be transformed into a system of 
ordinary differential equations (4.10) along characteristic manifolds, i.e. waves can be ob- 
served, travelling in three- and four-dimensional space-time domain along characteristic 
surfaces (a = 2) and volumina (a = 3), respectively. In Ex. 4.1, Ex. 4.4 and Ex. 4.7 (lin- 
ear case) and Ex. 4.2 and Ex. 4.5 (quasi-linear case) the wave propagation of mechanical 
systems, that have been introduced in Chap.2, will be analyzed. 


Example 4.1 (Wave propagation for the linear elastic bar formulation). The linear-elastic 
bar in one dimension (d = 1) has been introduced in Ex. 2.5. Following Sec. 4.1, the partial 
differential equation (2.14) can be replaced by the system of first-order partial differential 
equations (4.6) based on the state vector 


along with 
_ [pA 0 _ fo 4 fo 
paf ee-[2 3] reff oa 
The characteristics are then obtained as 
ds 
—=k'(t)=+c. 4.15 
=) = 40 (4.15) 


Accordingly, the slope of the characteristic lines in the space-time domain Q is determined 
by the constant velocity of wave propagation, already introduced in (2.16). In addition to 
condition (4.10), problem (4.10), gives rise to two further conditions 


d A , dx 
u 1 |_ k'(t)pA pATO| _ 
det | er k(t) 0 and det EA er 0. (4.16) 


Here, ö(t) = v(k(t),t) and n(t) = n(k(t),t) stand for, respectively, the velocity and the 
normal force along the characteristic lines s = k(t) whose slope k’(t) is given by (4.15). The 
two conditions in (4.16) eventually yield one independent equation of the form 
dô di 
k’(t)pA— - — =0. 4.17 
(HPA 7 (4.17) 
The system of ordinary differential equations (4.17) can be solved along the characteristic 
lines specified by (4.15). Alltogether, application of the method of characteristics converts 
the scalar second-order partial differential equation (2.14) into a system of four ordinary 
differential equations given by (4.15) and (4.17). These ordinary differential equations are 
supplemented with boundary conditions 


n(0,t) = P(t), o(L,t)=y'(), n(L,t) =0 (4.18) 


and initial conditions v(s, 0) = vọ and n(s,0) = 0. For the linear problem at hand, Eq. (4.15) 
and Eq. (4.17) are decoupled. Accordingly, Eq. (4.15) can be used in a first step to set up 
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the characteristic net consisting of straight characteristic lines specified by the constant velo- 
city of wave propagation c. A sample characteristic net thus obtained is shown in Ex. 4.23, 
Fig. 4.22. To each node (so,to) € Q of the characteristic net there is associated a nodal 
state vector zo = (vo,no). The unknown nodal states can be determined by integrating 
Eq. (4.17) along the characteristic lines. This procedure is illustrated with Fig. 4.1 in which 
three nodes Q and Pj, j € {1,2} are considered. The three nodes correspond to vertices on 
the characteristic net and thus comply with Eq. (4.15) in the sense that 

SB (-1)/e = 0 
io — Ip, 


is satisfied for j € {1,2}. The numerical integration of (4.17) relies on the finite difference 
scheme ne RER 

Der (ep a — = 
to — tp, to — Ip, 


J 


for j € {1,2}. Repeatedly applying this procedure to the characteristic net in Fig. 4.22 yields 
a (fully coupled) system of algebraic equations for the determination of the unknown states 
(va, na). Note that the states at the boundaries at s = L (cf. boundary conditions (4.18) 3) 
and t = 0 (initial conditions) are given. Eventually, boundary condition (4.18), determines 
the actuating forces f (tg) = ng for the nodes B lying on the boundary ats = 0. 


t Q 


Pi P, 


S 
> 


Figure 4.1.: Three nodes Q, Pı and P}, located on the characteristic net with associated states (vono), 
(up,,np,) and (vp,, np,). The characteristic lines kı and kz are associated with the slopes ds/dt = c and 
ds/dt = -c, respectively. 


Example 4.2 (Wave propagation for the geometrically exact string formulation). In line 
with the procedure outlined in Sec. 4.1, a curves = k(t) along that a solution zo (t) = z(k(t), t) 
is given, is called a characteristic curve if the partial derivatives 0,z and ðsz cannot be de- 
termined uniquely by the exclusive knowledge of z on curve k(t), i.e. the condition 


(E-K(t)D) az = F - DI) (4.19) 


A 0 0 I C 
p-[4 9] z-o] je (u) 
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gives rise to the condition 
det(E — k’(t)D) =0. (4.21) 


Consequently, k’(t) coincide with the eigenvalues of E relative to D (i.e. the eigenvalues 
of D'E). Correspondingly, there are 2d characteristics associated with the eigenvalues A; 
fori € {1,...,2d}. In particular, for the hyperelastic material model introduced in Ex. 2.1, 
function H in Eq. (4.20)2 assumes the form 


EA 1 EA 
H=—(|1-—|I+—p®p (4.22) 
2 v2 v4 


so that condition (4.21) yields the pairwise solution 


5 1 1 
k (t) = c 3 I + =) (4.23) 
ford =1, 
7 1 1 , 1 1 
k(t) = +c 3 1+ =} and k;(t) = +c 3 f = >) (4.24) 
ford = 2 and 


Al =) J1 1 
ki (t) = +c Al 1+ =} and k(t) = +c 3 I = =) (4.25) 


ford = 3. In Eq. (4.23) - Eq. (4.25) c denotes the constant velocity of wave propagation 
introduced in (2.16). 


Apart from the one-dimensional case (d = 1), the eigenvalues (4.24) and (4.25) are not all 
real-valued anymore for v < 1. This conforms with the fact that strings can sustain tensile 
forces but not compressive ones. The problem is then no more hyperbolic. 


Example 4.3 (Plane wave propagation for the geometrically exact string formulation). 
For the nonlinear string in two dimensions (d = 2), the second-order partial differential 
equation (2.11), is replaced by the system of first-order partial differential equations (4.6) 
based on the matrices in (4.20) and (4.22). The state vector z : Q > R‘is comprised of the 
vectorsv : Q +> R? andn : Q > R?. Then, two pairs of characteristics can be distinguished 
which are characterized by the slopes in (4.24). Correspondingly, two (i = 1,2) pairs of 
characteristics are defined by 


d z 
(=) = (-1)/*1¢; (4.26) 
dt} ;; 
where, with regard to (4.24) 
1 1 1 1 
cy =c 3 (i + >) and C2 =C 3 (i = >) è (4.27) 
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Each pair of characteristics has j € {1,2} members, where j = 1 refers to the forward 
propagating direction and j = 2 to the backward propagating direction. In addition to the 
directionality condition (4.26), the compatibility condition indicates a system of ordinary 
differential equations along the characteristics i, j € {1,2} given by 


dn ; dv 
U;-|— -1)/c;V; -|— -—1)/c;W =0 4.28 
Eee ee 
In Eq. (4.28), the functions 9 = v(k(t), t) andn = n(k(t), t) denote, respectively, the velocity 
and the contact force along the respective characteristic curve. Moreover, 


U;,=V;=P;p and W = pipe 


where 


P; = 64 i 0 + dj2 $ a X 


pı p 0 
Equations (4.26) and (4.28) constitute a system of coupled non-linear ordinary differential 
equations which can be discretized as follows: 


S P; R 
4 (-1Vei|p = 0, 
S (4.29) 
u;| TO ey + (-1)/(c:Vi)| 2O CR E (-1)/(ciW)| =0. l 
PU to -= tp; PU to ip, me 


In analogy to Ex. 4.1, vertices of the characteristic net are denoted by (so, tg), while the 
corresponding nodal states are denoted by (vg, ng). In contrast to Ex. 4.1, the characteristic 
net can not be determined beforehand, because c; in (4.29), depends on the state, see (4.27). 
Consequently, the location of the nodal net points, (sg, tg) and (sp,,, tp,,), in general need to 
be determined through the solution procedure, as well as the nodal states. 


ta Q 
ku kız 


Figure 4.2.: Integration along the characteristics k;; for i, j € {1,2} according to (4.29). The characteristics 
kij are associated with the slopes (ds/dt);; in (4.26). 


Scheme (4.29) is further illustrated with Fig. 4.2, in which the points Q and P;; located on 
the respective characteristics are displayed. This procedure can be repeatedly applied within 
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the space-time domain of interest. After assembly, a coupled system of non-linear algebraic 
equations is obtained which can be solved iteratively by means of Newton’s method. Prior to 
the solution, the boundary and initial conditions emanating from (2.11)23 need to be taken 
into account. Rem. 4.8 contains an outline of the implementation based on scheme (4.29). 


Example 4.4 (Wave propagation for the linear Timoshenko-Ehrenfest beam formula- 
tion). Consulting Rem. 2.14, the motion of a beam, formulated in terms of the Timoshenko- 
Ehrenfest formulation is governed by a system of second-order partial differential equations 
given by 


pAdu, -9,(H)-n=0, 
pla?@ — a,(M)+H-m=0. 
Regarding bending and shear deformations, constitutive relations are assumed as 
M=EI0,0=EIk ; H= GA(ðsu, + ©) = GAy. 


The initial boundary value problem at hand aligns with the framework proposed in Sec. 1.1 
forx = [u oe" along with 


(4.30) 


a-[2 al ae 0 


0 pl 0 4 ae eas 


ñ+ GAð; © 
H-m j 


Evaluating the directionality condition (4.10), gives rise to the characteristic directions 


Herein, cı and cz can be identified as the velocities of bending and shear waves, respectively. 
In Fig 4.3 the wave propagation for the Timoshenko-Ehrenfest beam formulation is illus- 
trated. Suppose, that Pıı Pız2 is an initial curve, along that the solution can be considered as 
given. Then the solution propagates along characteristic curves ki; for i, j € {1,2} such that 
the solution in point Q depends on points Pi; for i, j € {1,2}. 


ta Q 
kıı kı2 
Pir Pı2 
S 
> 


Figure 4.3.: Wave propagation for the linear Timoshenko-Ehrenfest beam formulation. 
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Example 4.5 (Plane wave propagation for the geometrically exact beam formulation - 
componential representation). Regarding Rem. 2.12 and assuming hyperelastic material 
(cf. Ex. 2.10), the function H : Q > R@ can be identified as 


ZA (1+yz?)cos® £4(1+y;*)sin© 0 
Hik = Bik + Op, Bijpj = -GAsin® GAcos® 0|. (4.31) 
0 0 EI 


where, use of 
EAy;*cos® for i=k=1 
Op, Bijpj = \EAy3’sin® for i=1Ak=2 (4.32) 
0 otherwise 


has been made. Note that, according to Ex. 2.8, the strain variable y is given by the relation 
Ya = ðsrı cos O + dry sin O. By evaluating the directionality condition (4.10),, the following 
characteristic slopes can be identified: 


EI\? GA\? 1EA(. 1\\? 
cı==+|1—| , &=+|— and c3 =+|-— {1+ = (4.33) 
pI pA 2 pA Y3 


Herein, c1, c2 and c3 can be interpreted as the velocity of wave propagation corresponding to 
bending, shear and compression, respectively. 


Figure 4.4.: Wave propagation for the Cosserat beam formulation for nonlinear constitutive relations. 


In addition to that, the compatibility condition indicates further ordinary differential equa- 
tions, governing the solution along characteristic lines. Together, a system of ordinary differ- 
ential equations is obtained that is equivalent to the underlying partial differential equation. 
This enables an integration of the problem at hand globally in space-time domain Q (cf. [31]). 
Observing that dsr = |losr|| [cos a sin a ol’, the strain variable ya might be reshaped 
as 

ya = ðsrı cos © + Acr2 sin O = (cos æ cos O + sina sin O) ||ösr|| 


= |lösr|| cos (a — ©) = |lösr|| cos £ . 
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Therein, use of the identities 


cos(@-®)+cos(@+®)=2cosa@cos®, 


cos(@-®)-cos(@+®)=2sinasin® 


has been made. See Fig. 4.5 for a geometrical interpretation. In case of assuming B = 0 
the velocity of longitudinal wave propagation of the geometrically exact beam formulation 
coincides with the velocity of longitudinal wave propagation of the geometrically exact string 
formulation (cf. Eq. (4.27) in Ex. 4.2). 


ðr 
dı 


D, 


D3 


Figure 4.5.: Planar illustration of the orientation of a cross-section of a geometrical exact beam. 


See Fig. 4.6 and Fig. 4.7 for an graphical evaluation of the longitudinal velocity c3 for 
the geometrically exact beam formulation depending on the shear angle p for different 
llösr|| € {0.5, 1.1, 1.4, 2.0, 5.0} and depending on |lösr|| for B € {0, 5,, 5. en, $n, T}, respec- 
tively. 


27 | 
1.7 at J 
el Be Er 


——— 5,0 

De tn 
Fi | 
er 
pee 

0.3 + 14 
[ | 
1 1 1 f 
0 n/4 n/2 32/4 T 
B 


Figure 4.6.: Longitudinal wave propagation c3 for the geometrically exact beam formulation, depending on £. 
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0.3 + 


Ilösr | 


Figure 4.7.: Longitudinal wave propagation c3 for the geometrically exact beam formulation, depending on 
llasrl- 


Example 4.6 (Objectivity). The velocities of the plane wave propagation for the geomet- 
rically exact beam formulation, which have been evaluated in Ex. 4.5, of course need to 
be independent to any coordinate transformation. Rotating the frame of reference, e.g. by 
A € SO(3), must not affect the characteristic directions. Hence, changing from componential 
to vectorial representation must lead to the same velocities of wave propagation. Thus, re- 
garding the formulation in vectorial representation and thereby taking the coefficients from 
Ex. 2.11, the matrix H can be determined as 


EA cos? © cosOsinO 0 
Hix = Bik + O,Bijpj = rat +yz”) |cos@sin® sin? © 0 (4.34) 
0 0 0 
sin? © -cos®sin28 0 0 0 0 
+ GA |- cos O sin 20 cos? © 0|+EI|0 0 0 
0 0 0 0 0 1 
In Eq. (4.34), the components of grad, (B) p are obtained as 
EAy;* cos’ © for i=k=1 
EAyz” sin? © for i=2andk=2 
2p, BijPj = (4.35) 
EAy;“cos®OsinO for i=landk=2Vi=2andk=1 
0 otherwise. 


where use of y3 = ðsrı cosO + Asr2 sin © has been made (cf. Ex. 2.8). By evaluating the 
directionality condition, the velocities of wave propagation are obtained as 


EI\? GA\? 1EA 1\\? 
cq =t+]—] . œ= +|— and c3 = +|-— |1+ — í (4.36) 
pI pA 
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Indeed a comparison of (4.36) with the velocities of wave propagation for the componential 
representation evaluated in Ex. 4.5 proves the objectivity. 


Example 4.7 (Plane wave propagation for the linear elastic continuum formulation). 
Subsequently we confine our attention to linear-elasticity. Assuming infinitesimal deforma- 
tions, the kinematical considerations that have been elaborated in Sec. 2.3 can be linearized 
in a direction u(s,t) around the reference configuration. Consequently, the geometrical lin- 
earization of the Green-Lagrange strain tensor is obtained as 


DaE(F(s + au(s, t)))|a=0 = ; (osu + (asu)") =E (4.37) 


Furthermore, assuming linear isotropy, the Cauchy stress tensor field ø (r,t) can be set con- 
stitutively into relation to the linearized strain tensor field e(r, t) as given by 


o = A(I : €)I + 2pe (4.38) 
giving rise to 
o =A(I: ðsw)I + u(ðsu + (Asu)!). (4.39) 


In Eq. (4.38) the two Lamé parameters A and p have been introduced. Inserting the divergence 
of (4.39) with respect to s, i.e. 


div ø = div(A(I : ðsu)I) + div(u(9su + (d,u)") (4.40) 
into the linearized balance of linear momentum 
divo+f=pou, 


a partial differential equation that governs the ininitesimal displacements u(s, t) of all ma- 
terial points s at any time t is obtained as 


(A + ðs (su) + uAu + f = pdru. (4.41) 


These differential equations are frequently referred to as ‘Lamés differential equations’ or 
‘Naviers differential equations’. Subsequently, these equations will be analyzed regarding 
the propagation of waves. Therefore, we restrict ourself to the planar case of linear elasticity, 
i.e. plane stress and plane strain is assumed, respectively. For more details on wave phenom- 
ena of solids we refer to [2] and [61]. For the sake of clarity componential notation is used 
subsequently. 


(i) Plane stress. Considering the special case of plane stress by assuming ciz = 03; = 0 
fori € {1,2,3} and £j3 = &; = 0 for j € {1,2}, the components of the Cauchy stress 
tensor are obtained as 

Oij = ÀEkk + 2uEij 


Consequently, a relation for £33 depending on £11 and £22 is given by 


£33:= (£11 + €22) « 


A 
2utA 
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This in turn leads to a relation for the components ofthe Cauchy stress tensor depend- 
ing on the linear strain components 


cij =A - Jekk + 2uei; VY i, j,k € {1,2} 


2u +4 


Consequently the divergence of the Cauchy stress tensor field yields 


A 
Oj Oi; safı _ sa 


O1(€11 + £22) +2 01£11 + O2€12 
02(€11 + £22) O1€21 + O2€22| ` 


Introducing the new functions p = [en E22 yel and v = [pu ru] and 
taking advantage of the property of mixed partials, i.e. 


OrE11 — 95,01 = 0, HE — As,02 =0 and Ayı2z - (9,02 + As,01) = 0, 


Naviers equations (4.41) can be established as a system of first-order partial differen- 
tial equations 
Ad, z + Bds,z + Cd;z = 0 (4.42) 


governing the solution z = [p ol". In Eq. (4.42) the coefficients 


a+2p a 0 0 0 0 0 H 0 0 
0 Ow 0 0 aat+20 0 0 
A=] 0 0 0 0 -1|, B=|0 0 0 -1 Of, 
0 00 -1 0 0 0 0 0 0 
0 00 0 0 0 0 0 0 - 

(4.43) 
0 0 0 -p 0 
a 0 0 0 -p 
C=|0 0 1 0 0 
100 0 0 
0 10 0 -1 


along with the constant parameter a = A(1 — A(2u. + 1)71), have been introduced. 
Evaluating the directionality condition 


det (C = ADız = BD;z) =0 
yields two distinct velocities of wave propagation 


2_ 2, 2_H Ze n ney eee 
cņ=c¢]1+c3=— and cp=ci+cn= 


1 4y(A +u) 


4.44 
p A+2p es 


Note that cz coincides with the longitudinal velocity of waves in an infinite plate 


(cf. [61, pp.79-83]). 
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(ii) Plane strain case. For the plane strain case, i.e. £i3 = &; = 0, Naviers equations (4.41) 
boil down to 


A E11 + dıE2 ci O2Y12 OF uy 0 
A +2 + -p| |= ; 
bee + oa F Es] F Eo P ur 0 
The problem can be written as a system of first-order partial differential equations 


taking the form of (4.42) along with the coefficients introduced in (4.43) by assuming 
a = À. The slopes of the characteristics can then be observed to be 


A+2 
acetal and cĉ =f +c = E. 
P P 
Herein, cr and cy, are representing the velocity of transversal and longitudinal wave 


propagation’, respectively (cf. [2, pp. 123-124] [43, pp. 318-320] and [61, pp. 10-15]). 


In Fig. 4.8 an illustration of the characteristic directions, along that longitudinal and transver- 
sal waves propagate in terms of planar linear elastodynamics of continua, are depicted for 
the planar formulation of linear elasticity of continua. 


£ a DD 
S2 82 
CT 
Sy 


Sy 


Figure 4.8.: Longitudinal and transversal wave propagation for the linear elastic continuum formulation - the 
plane stress and plane strain case. 


Remark 4.8 (Massaus scheme). In Sec. 4.1, it has been demonstrated, that the knowl- 
edge of characteristic manifolds enables a transformation of hyperbolic partial differential 
equations into a system of ordinary differential equations. Although ordinary differential 
equations are much more comfortable to solve than partial differential equations, complex 
mechanical systems lead to equations that are still challenging, even impossible, to solve ana- 
lytically. Therefore, a graphically-numerically approach, that could have been implemented, 
is introduced subsequently. 


Suppose, curve Jo is an initial curve, on that the solution z : Q > R” is considered as given. 
Then the solution at any point Q on curve J, depends on the solution at points P;; and can 


Waves that are traveling with velocity cz do not involve any rotation and waves that are traveling with 
velocity cr do not involve any dilatation. Thus cz and cr is frequently referred to as irrotational and equiv- 
oluminal velocity of wave propagation, respectively. Alternatively they are referred to as dilatational and 
distortional waves. 
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be determined by solving the system of ordinary differential equations (4.10) along charac- 
teristic lines kj; fori € {1,--- ,d} and j € {1,2}, by applying appropriate finite difference 
schemes. See Fig. 4.9 for an illustration thereof. Note that curves J, forn € {0,1,---} need to 
be strictly non-characteristic. Once the solution at point Q on curve J, has been determined, 
the solution among major characteristics, i.e. i = 1, can be interpolated. In the case of uni- 
directional wave propagation, interpolation is obviously not required (cf. e.g. Ex. 4.1). The 
given boundary and initial conditions can be applied directly at the corresponding nodes of 
the characteristic net. This graphical-numerical approach can be traced back to J. Massau 


(cf. [72]. 


Figure 4.9.: Illustration of the implemented graphical-numerical procedure - Massaus method. 


Remark 4.9 (Delay and differential flatness). Subsequently, some remarks on the delay of 
perturbations and its connection to the property of differential flatness are made. Therfore, 
the domain Q that is bounded by three non-characteristic curves Jo, Ko and K; as well as by 
the characteristic curve C} is considered. See Fig. 4.10 for a graphical illustration thereof. 


Delay. Assume the solution on curve K,(&) for E € [P2, Ps] is perturbed. Then, waves 
that are travelling along characteristic lines C will affect the solution on curve Ko(é) for 
€ € [Po, P5]. Since the solution on Ko(&) for E € [Po, Ps] in turn affects the solution on 
curve K;(&) for E € [Po, P4], the solution on curve K2(é) for E € [P3, P4] is affected by 
reflecting waves caused by perturbations imposed on curve Kz (E) for & € [P2, P3]. Note that 
the solution on K2(&) for ë € [Pı,P2] is governed by initial conditions that are considered 
as given on curve Jo. In the context of the inverse dynamics problem, we refer to curve 
K2(&) for ë € [Pi, P2] as pre-actuation phase and to curve K2(é) for € € [P3, P4] as post- 
actuation phase. 


The span of both pre- and post-actuation phases are determined by the slope of the char- 
acteristic curves together with the distance between the boundaries. The shaded domain in 
Fig. 4.10 elucidates this correlation. (See also Ex. 4.23, Fig. 4.22 for the corresponding linear 
problem.) 
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Figure 4.10.: Characteristic net. 


Flatness. Concentrating on the space-time tube between the two characteristic curves Cy and 
Cj the solution z(Kı) at all points on curve K, depend on the solution z(K2) at points on K3. 
Once the solution z(Kı) at all points on K, has been determined, the solution z(Ko) on Ko 
can be established by repeating this procedure. Thus, the solution z(Ko) on Ko can be totally 
parameterized by the solution z(K2) given on K3. 


This is in analogy to the property of differential flatness introduced in Def. 3.23 and accom- 
panying Ex. 3.24, Ex. 3.25 and Ex. 3.26. 


Remark 4.10 (Riemann Invariants). In general, the solution ë € R +> z(£) itself is not 
conserved along characteristic lines. But there might exist quantities, that are preserved 
along these lines. This remark aims to reveal such invariant functions. 
A function z € R°4 + h(z) € R@ is a preserved quantity along characteristic lines if 

öch(z(2)) = 0 (4.45) 
holds. Together with equation (4.6) 

ðgz = Qz Det + OZ Des 

and setting therein for simplicity F = 0 the following relation can be established: 


(-(D" HT. + ciðzh) ðz=0. 


Thus, the eigenvectors ðzh associated to the eigenvalues of D’'E yield along characteristic 
lines the preserved quantity h(z). The following eigenvalue problem can then be established: 


(D'E) əzh = ciðzh . 
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Finding the eigenvalues c; and the corresponding eigenfunctions, the preserved function h(z) 
can be found by integrating the eigenfunctions. The function h(z) is called Riemann Invari- 
ant of the wave equation at hand (cf. [29] or [7] for more details). For the linear problem, 
introduced in Ex. 2.5 and analyzed in Ex. 4.1, the Riemann Invariants can be determined as 


follows. 


; 1 
The eigenvectors p; associated with the eigenvalues c; = (—1)'*! (Ep~')? of the matrix 


=1p\T _ 0 -EA 
MET Poe 0 | 
can be obtained as 5 
pi = (%h);= |-pAc; 1] . (4.46) 


Integrating these eigenvectors (4.46) yields the Riemann Invariants 
h;=n- pAcjq = const. (4.47) 
See Ex. 4.23, depicting both preserved functions of the problem at hand. 


Remark 4.11 (Hodograph Transformation). Alternatively to Rem. 4.10, Riemann Invari- 
ants might be found by applying a hodograph transformation to the partial differential equa- 
tion at hand and searching subsequently for characteristic lines. The hodograph transforma- 
tions is given by F(z) : RE + R@ such that x > F\(z) (cf. eg. [29]). For instance, 
system (4.6) with (4.14), can be transformed by inserting the Jacobian 


OgS Ons 
Ozx =|” 
A 2 T 


along with its inverse 


o%q q 


= (0x) = (OgSOnt — Onsögt) | Int | 


—Ogt qs 
into the system of first-order partial differential equations 
D(z)dgx + È(z)ðnx = F(z) . (4.48) 


In Eq. (4.48) the following coefficients have been introduced: 


a fo a]l Ss: SpA. 0 Jo 
b=] sepet 8) eft. a 
Evaluating the directionality condition leads to the desired Riemann Invariants. A line 
q = k(n) is a characteristic line if 


[e-5%,) 
det | E - D— | = 0, 
dq 
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hence 
2 E 
(Dgn)* = pAgEA e Dgn = +pAo4|— . (4.50) 
p 
By integrating Eq. (4.50), the Riemann Invariants are obtained as 
h;=n-cipAog = const. (4.51) 
Compare Eq. (4.51) with Eq. (4.47). 


Remark 4.12 (Geometrical interpretation: The general non-linear case). To gain more 
insights into hyperbolic problems in general, it might be beneficial to reflect some geometrical 
aspects of general first-order non-linear partial differential equations as given by 


F(s,t,z,p,q) =0. (4.52) 


Therein p = 0;z and q = 9ız have been introduced denoting the partial derivatives of the 
solution z(s,t): Qh Z C R with respect tos e S C Randt E T C R whereQ=SXT. 
Solving Eq. (4.52) basically implies to find an integral surface z(s, t) that satisfies Eq. (4.52) 
along with given boundary and initial conditions. In analogy to the basic concept of direction 
fields within the theory of ordinary differential equations (cf. [12] and [92]), the partial 
differential equation (4.52) allocates to every point P(sp,tp,zp) € A= SXTXZ c R 
planes 

(sp —s)p + (tp — t)q — (zp - z) = 0, (4.53) 


that are tangential to integral surfaces of Eq. (4.52). Thus, Eq. (4.52) inheres the orienta- 
tion of the tangent planes of the solution z(s,t) at each point P € A. Since we claim that 
(dpF)* + (dgF)* # 0, Eq. (4.52) forms a one-parametric family for the partial derivatives 
Am p(A) and à q(A) such that the planes (4.53) are envelopping surfaces of cones (see 
Fig. 4.11). To each point P € A one cone is allocated, touching the integral surface. In the 
special case of (quasi-) linear partial differential equations these cones degenerate to single 
lines and the envelopping surfaces form a pencil of planes (cf. Rem. 4.13). 


w(A) 


P(sp, tp, zp) 


t s 


Figure 4.11.: Integral surfaces envelopping Monge’s cone. 
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These cones are referred to as Monge’s cones, and its generating lines are spanning a one- 
parametric vector fieldw : Rt A as given by 


pF 
w(A) = OgF 
PoOpF + qdgF 
It is called Monge’s vector field and curves k(£) whose tangents Dek(&) belong to Monge’s 


vector field, are referred to as focal curves or Monge’s curves. This may be achieved by 
comparing the derivative of Eq. (4.52) with respect to A 


DaF (s, t, z, p(A), Q(A)) = OF Dap(A) +94 FDig(A) = 0 


as well as the derivative of the envelopping surfaces, introduced in (4.53), both with respect 
to À 


(sp - s) Dap (å) + (tp - t) Dig(A) = 0 
such that the following relation of the plane increments can be established: 
(sp — s) : (tp — t) : (zp — z) = OF : OGF : (pdpF + qdgF). 


Ifz(s, t) is a solution of (4.52), then Monge’s curves are obviously part of the integral surface. 
Consequently, they are called characteristic curves. 


Figure 4.12.: Monge’s cone with integral surface of initial boundary value problem along characteristic curve 


Kid). 


Along these characteristic lines, the partial differential equation (4.52) can then equivalently 
be replaced by a system of ordinary differential equations as given by 
s Op F 
Dek(&) = De] t| = OgF =w. (4.54) 
z DPF + qogF 
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Alternatively, sytem (4.54) might be reformulated as 


Ds=— Diz =q 2 p (4.55) 
s=—, z=gqtr—p)D. 5 
i ; OqF 


In addition to that, the set of all tangent planes along Monge’s curves are referred to as 
Monge’s strips. In case of z(s,t) is a solution to the differential equation (4.52) then these 
strips are called characteristic strips. These characteristic strips give rise to two further or- 
dinary differential equations revealing the change of the partial derivatives q and p with 
respect to č 


Dep = —-0;F — pdzF, Deq = -9,F — qozF (4.56) 
or alternatively with respect to t 


IF  a,F OF aF 
ER Dıq = -— - 
qF 


(4.57) 


This is obtained by analyzing the derivative of Eq. (4.52) with respect to s 
DsF(s, t, z(s, t), p(s, t), q(s, t)) = 9sF + pdzF + 3pFðsp + 0,Fd,q = 0. 


By taking additionally advantage of the equality of mixed partials, i.e. 0%,z = 3% z eg. 
ðsq = O:p, such that the following relation can be established 


OF + p0,F + OpFOsp + 3gFðp = 0. (4.58) 


Then Monge’s vector field (4.54) gives rise to dyF = Dzs and 0,F = Det, hence Eq. (4.58) 
assumes the form 


O;F + po,F + Desðsp + Deto;p = ðsF + pozF +Dep =0. 


Herein use of Dep = (Des) (Asp) + (Det) (drp) has been made. Accordingly, by taking the 
derivative of Eq. (4.52) with respect to the variable t, the following holds: 


D;F(s, t, 2(s, t), p(s, t), q(s, t)) = 0;F + q0,F + Desd,q + Detd;q = OF + qo,F + Deq =0. 


For more details on the geometry of partial differential equations, see [84, pp.41-49] , [31], 
[11, pp.369-370], [33] or [37]. 


Remark 4.13 (Geometrical interpretation: The special quasi-linear case). Based on the 
results proposed in Rem. 4.12, the special case of quasi-linear partial differential equations, 
assuming the form 


F(s, t, z, 952, dz) = a(s, t, z)0;z + D(s, t, z)ðsz — c(s, t,z) =0, (4.59) 


is analyzed subsequently. Quasi-linearity implies linearity in the highest occuring partial 
derivatives, i.e. the coefficients a,b and c in (4.59) are allowed to depend on the space and 


86 


4.1. Wave propagation 


time variable as well as on the solution itself. Consequently, the coefficients are forming a 


vector field 
w(s,t,z)=[a b cl’ (4.60) 


that is perpendicular to the (normal) vector field 
n=[az az -1]'. 


Hence, w(s, t, z) is tangential to the solution z(s, t) in every pintPe A=SxTxZcC R3. 
The vector fieldw(s,t,z) is the Monge’s vector field, introduced in Rem. 4.12. In each point 
P € A Monge’s vector field forms an axis of intersecting planes that are tangential to the 
integral surfaces. The orientation of Monge’s vector field gives rise to a system of ordinary 
differential equations given by 


ds:dt:du=a:b:c. (4.61) 


These equations are equivalently governing the solution of (4.59) along characteristic curves. 


P(sp, tp, zp) 


Figure 4.13.: Integral surface (contact element) z — Zo = po(s — so) + go(t — io) at point Po with normal and 
tangential vector field n and w, respectively. 


Considering an initial line along which the solution is given and is everywhere tangential to 
Monge’s vector field w it becomes evident, that data which is considered as given only along 
such lines does not suffice to obtain an unique solution for the tangential planes. Hence, 
discontinuities of the partial derivatives d,z and 0;z, propagating along these characteristic 
lines, need to be permitted strictly. Consequently, 


dz = 0,z ds + 9,2 dt ie. 02 = z" — d;zk'(s) (4.62) 


along with the given quasi-linear partial differential equation (4.59) lead to an indifferent 
solution for the partial derivatives ðsz and dız. 
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Figure 4.14.: Monge’s pencil of planes. 


Example 4.14 (Linear transport equation). In the simplest case, linear transport phenom- 
ena are governed by partial differential equations in form of 


F(p,q) =ap+bq-c=0, (4.63) 


where a, b and c are assumed to be constant. The partial derivatives of the solution z with 
respect to both independent variables s and t are defined as p = 0,z andq = 02, respectively. 
In accordance to Rem. 4.12 and Rem. 4.13, the differential equation (4.63) gives rise to a 
parameterization of the partial derivatives p and q, as given by 


q=<cA and p=<(1-A), (4.64) 
b a 
such that Monge’s vector field 
Op F a 
w= qF = |b 
PoOpF + qogF c 
can be identified as the axis of a pencil of planes given by 


z- zp = (s — sp)p + (t — tp)q = (s SEE zit ZA: (4.65) 


Thus, along characteristic lines, the solution is governed by a system of ordinary differential 
equations given by 


D ee d D ee + C 
s= — =- an a —p= —p=-. 
DE TAS ap 


Due to the linear structure of the underlying equation, the partial derivatives q and p are 
only changing linearly with respect to s and t, i.e. geometrically the characteristic stripes 
are not distorted. Consequently, 


D;p=0 and Diq =0. 
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Example 4.15 (Eikonal equation (cf. [37])). Maybe the simplest representative of general 
non-linear first-order partial differential equations is the Eikonal equation given by 

F(s,t,z,p,q) =p +q -1=0. (4.66) 


According to Rem. 4.12 and Rem. 4.13 the corresponding Monge’s vecor field takes the form 


OpF p p 
w= OqF =2 q =2|q 
PoOpF + qogF p + q? 1 


By introducing the parameterization 
p=cosA and q=sind 
a one-parametric vector field w(A) : R > A can be obtained as 
cos A 
w(A) =2]sinda (4.67) 
1 


determining generating lines for Monge’s cones. This cones are envelopped by the tangential 
planes 


(z -= zp) = p(s — sp) + q(t — tp) = cosA(s - sp) + sin A(t — tp) . 


0,5 —0.5 


Figure 4.15.: Normalized Monge’s vector field for the Eikonal equation (4.66). 


Along characteristic lines Eq. (4.66) can equivalently be transformed into a system of ordi- 
nary differential equations given by 


pF 2 : 
Dis = — = = and De le (2) }. 
q 
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Although, the Eikonal equation (4.66) is a non-linear equation, Monge’s cones do not change 
its shape in the domain A. Consequently, Monge’s strips are not distorted, i.e. 


Dip =0 and Diq=0. 


Example 4.16 (Hamilton-Jacobi equations). Another prominent representative within the 
class of general non-linear partial differential equations is given by 


F(r,t,S, p,q) = 4S + H(r,t,9,5) =q+H(r,t,p) =0. (4.68) 


This non-linear first-order partial differential equation, governing the action 


s= fu, r,D;r) dt , 


is called the Hamilton-Jacobi equation. In Eq. (4.68) the function H denotes the Hamiltonian 
of the system and the partial derivatives with respect to time t € T and configuration r € 
R"d are defined as 

q=aS and p=as. 


Following the geometrical interpretation of general non-linear partial differential equations, 
analyzed in Rem. 4.12, the Hamilton-Jacobi equation (4.68) can be transformed into a system 
of ordinary differential equations along characteristic lines. Evaluating (4.54)/(4.55) yields 


OpF 
Dr = == =H, (4.69) 
OqF 
pF 
D:S = q+ p— =qg+p%H. (4.70) 
OgF 


Due to the non-linearity, the partial derivatives q and p are not changing linearly anymore. 
Hence, along characteristic strips, (4.56)/(4.57), leads additionally to a system of ordinary 
differential equations given by 


OF OsF 
Dip = - —- p— =-9,H, 4.71 
tP OqF Par r ( ) 
aF OsF 
Daq = -+ -q> = -3H . 4.72 
tq agF ER; t ( ) 


The two ordinary differential equations (4.69), and (4.71), can be identified as the well- 
known canonical Hamiltonian equations. These equations are the characteristic directions 
of the Hamilton-Jacobi equation. Note that the action S does not occur explicitly in (4.68), 
such that dsF = 0. Furthermore ö,F = 1 applies. For more details on the Hamilton-Jacobi 
theory, see e.g. [11, pp. 248-270], and [84, pp. 66-68]. 


Example. Assuming natural, autonomous mechanical systems with Hamiltonian as given 


by 
H=p°+V(r). (4.73) 
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Then, the Hamilton-Jacobi equation (4.68) is obtained as 

F(r,t,S,p,q) =q+p?+V(r)=0. 
In accordance to Rem. 4.12 and Rem. 4.13, Monge’s vector field is determined as 

Op F 2p 
w= OgF =| 1 l. (4.74) 
PopF + qogF 2p? +q 

By introducing the parameterisation 

par and q= -1? - V(r) 
Monge’s vector field constitutes a one-parametric vector field 

2A 


w(}) = 1 (4.75) 
A? -V(r) 


that gives rise to the generating lines of Monge’s cones envelopped by the tangential planes 
(S = So) = Ar = ro) - (4? + V(r))(t = to). 


Along characteristic lines the Hamilton-Jacobi equation (4.68) can then equivalently be 
transformed into a system of ordinary differential equations as given by 
D 2 d DS 2 2p? 
r=— = an =qt—p=qt : 

t agF P t =q a re a 
In contrast to Ex. 4.14 and Ex. 4.15 for the Hamilton-Jacobi equations Monge’s cones change 
its shape in any point A, consequently Monge’s strips are distorted. This leads to an addi- 
tional system of ordinary differential equations given by 

O,F oF 


D = = = -D,V d Dig = -— = 0 
p= =-DV) and Dg =-S 


Figure 4.16.: Normalized Monge’s vector field for the Hamilton-Jacobi equation (4.68) assuming the Hamilto- 
nian (4.73). 
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Remark 4.17 (Plasticity). The elasto-plastic behaviour of general continua is governed by 
quasi-linear hyperbolic partial differential equations as well (cf. e.g. [84, pp. 117-118] and 
[1]). Following [83] and [97], the characteristic lines can be identified coinciding with the 
directions of plastic flow. 


4.2. Space-time finite element method? 


In Sec. 4.1 the hyperbolic structure of the governing equations has been revealed and 
consequences regarding the mechanisms of how information flows in physical time-de- 
pendent systems have been elaborated. The method of characteristics, proposed in Sec. 
4.1 as well, relies on an integration of the partial differential equation simultaneously in 
space and time. This intuitively motivates the development of further integration meth- 
ods, that are based on a simultaneous discretization of the space-time domain. Therfore, 
the present section aims to present novel Galerkin space-time finite element methods 
for the inverse dynamics of infinite-dimensional systems. For further informations on 
space-time finite element methods in general we would like to refer to the two pioneer- 
ing publications [9] and [52]. See also [80], [6], [54], [53] and [47]. 


By introducing the velocity v(s, t) = 0;x(s, t), the underlying partial differential equation 
at hand (1.1) can be transformed into a system of partial differential equations that is 
first-order in time: 

a,x-v=0, 


4.76 
Ad;v — 0;(Bd,x) =C. (278) 


Multiplying the equations in (4.76) by sufficiently smooth test functions w1, w2 : Q +> R@ 
and subsequently integrating over the space-time domain Q = S x T yields 


[jr axo) dQ =0 
2 (4.77) 

f» - (Að;v — ðs (Bəðsx)) dQ = f» -C dQ 
Q Q 


Integrating by parts the second term on the left-hand side of (4.77), and taking into ac- 
count the boundary conditions (2.11), yields 


[ome Aa das f awa: Bax dQ + FE) wal(styeaa, dt 
Q Q 89 f 


= f» -C dQ +f n(t) E wal(sneoo, di. 
Q 


09 


(4.78) 


3 This section partly reproduces [93]. 
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Concerning the prescribed boundary condition on Q}, we make use of weakly enforced 
servo-constraints of the form 


f wi) (VO ~ slenca) at = 0 (4.79) 


90 


where w; : 9Q, > Rf is a third test function. 


Now the weak formulation of the inverse dynamics problem at hand can be stated as 
follows: Given b, fy, xo and vo, find x € V),v € V, and f € V3, such that 


[om Alax—0)d0=0 
Q 


[ow - Ad,v dQ +f Oswa : Bd,x dQ + fi): Wal (st) cane dt 
Q Q 89 i l 
(4.80) 


= [wcag f TETE 
Q 


ay 
f mO WO- kane) dt =0 
any 


for arbitrary test functions we € Wg for a € {1,2} and w3 € Wz. The prescribed initial 
conditions (2.30) give rise to the Dirichlet-type boundary 90, of the space-time domain 
Q. Accordingly, we have 


V = {x : Qt RI |x = xo on AQo} 

V = {v : Qh RI |v = vo on AQ} 

V = {f:3Qf m Rİ |f = fo on 3f N aN} 
and 

Wa = {wa : Q m RË |wa = 0 on aN} 

W = {ws : 90, > Rf |w; =0 on IQy N ANo} 


fora € {1, 2}. The space-time finite element discretization of the weak form (4.80) is based 
on finite-dimensional sub-spaces vi C V; and wh C W; associated with interpolations of 
the trial functions 


x'(s,t) = Y Nals, t)xa fit) = E LOS: 
a and h iENo 
v”(s,t) = x Nals, t)Va n) = 2, Ln; 
where no is the total number of nodes giving rise to the index set Mo = {1,...,no}. 


Moreover, index set Nf C No contains the nodes lying on dQ. Similarly, N, C No 
contains the nodes lying on 0Q,,. The test functions are approximated by 


no 
w! (s, t= X Nals, Waa and w(t) = > L;(t)ws, 
a=1 ieNy 
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for a € {1,2}. 


Remark 4.18. (Direct dynamics) The space-time finite element method has been introduced 
in Sec. 4.2 for the inverse dynamics problem. Therefore, the servo-constraint has been im- 
posed weakly. In principle, ideal-orthogonal (contact-)constraints can be imposed in the same 
manner, such that the proposed space-time finite element method is also capable of solving 
the direct dynamics of constrained mechanical systems. 


Resulting algebraic formulation. 


In the present work, we confine our attention to standard low-order isoparametric finite 
elements. In particular, L; denote a-linear Lagrangian shape functions also employed 
in Chap. 3, while N, denote («æ + 1)-linear Lagrangian shape functions. Inserting the 
above finite element approximations into weak form (4.80) yields the algebraic system of 
equations 


AQ- MV =0, 
AV + Fi" (Q) + BY Fr = Fxty B, Fy , (4.81) 
B,O -Tr=0 


where 


Aab = u f ANað:Np dQ , 
Q 

Ma = nal ANN, dQ, 
Q 

Fa = | a,NaBa;x" dO , 

Q 

Br, = u f LiN, dt, (4.82) 
OQ ¢ 

Bri = u f LiNp dt 3 
aQ 


PPE [mc dQ, 
Q 


r= f Liy dt. 
AG; 


In (4.81), the unknown nodal quantities xa, vg and f; have been assembled in correspond- 
ing nodal system vectors Q, V and Ff. Elimination of the nodal velocity vector V yields 
the residual vector 


AM AQ + F™(Q)+ B Fp- F™ -B Fn 


8.0-7 (4.83) 
n 


R(Q, Ff) = | 
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The resulting algebraic system of non-linear equations, R(Q, Ff) = 0, can be solved 
iteratively for the nodal unknowns (Q, F f) by applying Newton’s method. In this con- 
nection, the iteration matrix is given by 


AM 'A+DF™(Q) By 


B, f (4.84) 


DR(Q, Ff) = | 


where DF™ (Q) results from (4.82)3 together with the constitutive term (2.12)2. A straight 
forward calculation yields the nodal contribution to DF (Q) given by 


dx, F=" (Q) = i 0;NaHa;Np dQ (4.85) 
Q 


where matrix H has been defined in (4.5). Note that p therein has now to be substituted 
with x". 


Recursive implementation. 


The newly devised space-time finite element method yields the system of algebraic equa- 
tions (4.81) which can be solved by applying the iterative procedure outlined in Sec. 4.2. 
This implies the simultaneous solution of all the unknowns resulting from the space-time 
discretization. 


T | rear ta ay SS ey 
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! i H H 

i | | l 

en bosi 

| | i ! 

H | I 1 
nor A | Mn 

I H H H 

I 1 z 1 

[i 1 I H 

Tü j p | 

Ka Ed 

H 1 i ! 

| H H H 
Dane Ge Se ee eee i = 
So=0 S Sn-1 Sn Sn-ı Sn=L S 


Figure 4.17.: Division of the space-time domain Q = S x T into time-space slabs Qn = ($n-1,5n) X T for 
n=1,...,N. 


We next show that the solution can alternatively be attained by applying a recursive 
scheme which relies on the decomposition of the space-time domain Q into N time-space 
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slabs Qn = (5-1, Sn) X T, n = 1,..., N (cf. Fig.4.17). Focusing on one representative slab 
Qn, the recursive solution procedure relies on the weak form 


i w,-A-(&%x—-—v) dQ=0 
Qn 


f wA aw doe | dwa: Bax d+ | w2: fn- dt 
Qn n In-ı 


= [ we Cdr f wa: fy at (4.86) 
Qn Tr 
fw e-r) dt = 0 


In 


J W3n-1 g (x ~ Visi) dt =o 
In-ı 


where Ip = {sn} X T. Note that Ty = 0Q, and Ip = 90; is defined. In essence, weak 
form (4.86) emanates from the original formulation (4.80) by restricting the space-time 
domain to Q,. In this connection, (4.86)3,4 have been introduced to facilitate a convenient 
description of the recursive solution procedure. The overall space-time approximation 
remains the same as before. The approximation of the newly introduced intermediate 
quantities f„ (n =1,...,N - 1) is given by 


TOS E LOLs, (4.87) 
iENn 


where index set Mn contains the nodes lying on Iņ. In this way, fi € ye C V3,, where 
Vs, = {fai Tn ORF, = Fn on Tn Qo} 

and f „(n=1,...,N - 1) have to satisfiy the initial conditions. Similarly, we introduce 
W3, = [ws, : Ip |> R? |w3, = 0 on Ip N a2) 


for n = 0,..., N — 1, so that in (4.86)34, w3, € W3„. Correspondingly, w € wi C W3, is 
given by 

wy = 3 L;(t)(w3,)i 

iENn 

It can be easily shown that the recursive application of (4.86) for n = 1,...,N is equiva- 
lent to the original space-time method, provided that the approximation of y, in (4.86)3,4 
conforms with the continuity of the displacement approximation. To ensure this, we 
choose 

yta) = > Li(t)y,, (4.88) 

ieN„ 
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so that (4.86)34 enforce rl; = yh and rh 35 yha Now, the equivalence between 
(4.86); and (4.80); follows from the property 


ARE da= fe) dQ. 


Similarly, the equivalence between Eq. (4.86) and Eq. (4.80), follows by additionally 
taking into account 


N 


D| weh ar- [ wef, at) = [min ar- [weh de. 


n=1 


Moreover, (4.86); contains (4.80); for n = N under the provision that the prescribed tra- 
jectory y on Tzr is based on nodal interpolation of the type (4.88). 


Now the proposed space-stepping algorithm results from the recursive application of the 
discretized weak form (4.86) in backward space direction, that is, for n = N,N — 1,...,1 
as shown in Tab. 4.1. 


LOAD y}, fh 
DO n=N,--- ‚1 


SOLVE Eq.(4.86) > FR» oo» vi ft 


n-1 


Table 4.1.: Recursive implementation. 


Note that for n = N the given data follows from the prescribed data on Tz leading to 
yn, = y" and fi = fe where the prescribed data y and f; is assumed to be given in 
interpolated form of the type (4.88) and (4.87), respectively. 


The implementation of the recursive scheme is again based on the algebraic formulation 
outlined in Sec. 4.2. However, in each step of the recursive scheme all of the algebraic 
quantities are now confined to the respective slab Q, (see also Rem. 4.20). 


Remark 4.19. The meaning of the newly introduced intermediate quantities f,, and f „_ı 
(n=2,...,N — 1) in Eq. (4.86)2 can be elucidated by considering the Euler-Lagrange equa- 
tions emanating from Eq. (4.86)2. A straightforward calculation based on integration by 
parts applied to the second term on the left-hand side of Eq. (4.86)2 yields 


J wa : (Adv — ðs (Bd,x) — C) da [ wa: (n-f,) ars f w2: (fa -n) dt = 0 
Qn 


In Th-1 


Accordingly, f, = nlp, and f„_ı = nlp,_,, so that f,, and f,„_, correspond to the contact 
forces in the string along I„ and Tn-1, respectively. 
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Remark 4.20 (Numerical effort). Assuming a regular discretization of the space-time do- 
main Q = S xT relying on ns nodes in space direction and n; nodes in time direction 
amounts to a total number of nodes ng = ns : nı. The simultaneous solution for all un- 
knowns outlined in Sec. 4.2 relies on the iteration matrix (4.84), whose dimension is equal to 
(ns+1)- (nı-1):d. In contrast to that, the dimension of the iteration matrix corresponding to 
each of the N = ns—1 steps of the recursive solution procedure is equal to (2a+1)-(n—1)-d. 


Remark 4.21. The semi-discrete approach described in Sec. 3.1 can be linked to the space- 
time finite element method developed above. To this end, we consider a Galerkin-based dis- 
cretization in time of the semi-discrete formulation. In particular, we introduce the following 
approximations of the nodal position and velocity vectors 


qi (t) =L;(t)g;; 


ob (t) = LHi; (4.89) 


For conciseness the summation convention applies throughout this remark. Moreover, for 
simplicity we again employ Lagrangian shape functions in (4.89). In addition to that, we 
introduce nodal weighting functions of the form 


Wi(t) = Ma(t)Wia (4.90) 


where My are basis functions to be specified below. Now the semi-discrete equations (3.51) 
can be recast in the weighted Galerkin form 


fw -Mij (Daqt -— vh) dt=0 


JWE- (Ma Deoh FA -FO + 580 -ipa fiO) dt=0 (491) 


fw ` Ôi p+ (a! - y) dt=0 


Note that (4.91); serves the purpose of introducing the nodal velocities vh. Time-stepping 
schemes typically applied in the context of the semi-discrete formulation can be recovered 
from (4.91) by choosing specific discontinuous shape functions Mı in (4.90), along with 
associated quadrature formulas for the evaluation of the time integrals. This approach gives 
rise to the so-called continuous Galerkin method! (see [36]). For example, focusing on the 
first two equations in (4.91) emanating from the underlying system of first-order ordinary 
differential equations, choosing constant shape functions Mı and linear Lagrangian shape 
functions L;, along with the mid-point quadrature for the evaluation of the time integrals 


yields the mid-point rule (cf. e.g. [19]). 


4 Note that the designation continuous Galerkin is attributed to continuous test functions despite the use of 
discontinuous weighting functions. This is in contrast to the so-called discontinuous Galerkin method which 
relies on test and weighting functions, both being discontinuous (cf. [36]). 
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In sharp contrast to that, the newly developed space-time finite element method can be re- 
covered from (4.91) by choosing continuous shape functions Ma. In particular, using again 
Lagrangian shape functions, i.e. Ma(t) = La(t), (4.91) leads to the algebraic formulation 


(4.81), 


(i) 


(ii) 


provided that 


the shape functions Na(s, t) used in the above space-time finite element method have 
tensor product form (cf. [16, 81], see also [99]). That is, 


Na(s, t) = Li(t)L; (s1)Lr (s2) - Lm(sa) - (4.92) 


Referring to the local node numbering of a a-linear master element the node numbers 
in (4.92) are related to each other. In Fig. 4.18 and Tab. 4.2, this relation is illustrated 
exemplarily for a = 1. 


4 ji 3 

1 2 
Table 4.2.: Relation between the local node num- Figure 4.18.: Bi-linear master element with local 
bers of the shape functions in (4.92). node numbering. 


In this connection, the one-dimensional linear Lagrangian shape functions are given 


by 


(4.93) 
12(8) = 


where & and & refer to the nodal coordinates, while € stands for eithers € Sort € T. 


rectangular shaped space-time finite elements are used leading to a structured mesh 
in space-time. Note that this precludes unstructured space-time meshes which are 
feasible with the space-time finite element method proposed above. 


To further illustrate the connection between (4.91) and (4.81), we consider the first term in 
(4.91)2, which yields 


[wh mde! a wu f La(t)L,(t) dt Mijo jk 
T T 
-Wa [rane dt fano) dsv jk 
s 


= Wia [Ap Lj(s)La(t) Lj(s)L,(t) dQ vx 


Wa Na(s,t) OrNp (s,t) Ve 
= W2, ` AabVb 
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Here, M;; introduced in (3.7) has been used. Moreover, the relation between the nodal shape 
functions in Tab. 4.2 has been taken into account. 


The remaining terms in (4.91) can be linked to the corresponding terms in (4.81) in a similar 
manner. 


Remark 4.22. [A discontinuous Galerkin approach] Due to the hyperbolic structure, the so- 
lution of the partial differential equation is dominated by wave phenomena. The mechanisms 
of wave propagation have lead to the perception that a simultaneous space-time integration 
might be the natural choice of solving the inverse dynamics problem of infinite-dimensional 
systems. Therefore a modified discontinuous Galerkin method is derived below. Its virtue lies 
in not being confined to spatially continuous systems. We are considering ordinary differen- 
tial equations in form of 

D;z+F(z,t)+Uu =0 (4.94) 


governing the solution z = lq o|” that is additionally subjected to servo-constraints as 
given by 
h(q,t) = Hq - y(t) =0. (4.95) 
In Eq. (4.94) 
0 
U = 4.96 
ar 
has been introduced. Sticking to classical literature on (discontinuous) Galerkin methods 
for ordinary differential equations (cf. e.g. [65], [35] and [36]), an approximate solution 
z" to (4.94) is searched for in a finite dimensional subspace spanned by polynomial basis 
functions pk (Ta) on subintervals Ta = [ta,tarı] C T. Herein, Pk denotes the space of all 
polynomials of degree < k. For the sake of clarity, a uniform discretization tay, = ta + Ah 
foro < A < N and Nh =T is assumed. Then the approximate solution z” is searched for, 
such that 


I w (Dız(t) + F(z,t) +U) dt +w4 [za] +wa,,- Hi zau] =0 (4.97) 
A 
is satisfied for allw € P*. In Eq. (4.97) the abbreviations 
aS m z(ta+č) and z, = u z(ta — &) (4.98) 
é-0 20 


has been introduced, denoting the limit ofz from above and below, respectively. Furthermore, 
the bracket [za] = z(t})-z(t) represents the jump ofz at time ta. An illustration thereof 
is depicted in Fig. 4.19 fork = 1. 
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ta t 


Figure 4.19.: Discontinuous approximation of the test and trial functions for k = 1. 


The spaces of the test and trial functions w and z are chosen to be the same. In general, 
test and trial functions are discontinuous at the nodes ta. Consequently, initial conditions 
z(t = 0) = Zo need to be imposed variationally by adding 


wa: [zal 


to Eq. (4.97). A strong imposition would lead to an overdeterminacy of the coefficients of z. 
Additionally, the servo-constraints (4.95) are imposed variationally 


f» . (Hq - y(t) dt=0 (4.99) 
T 


by adding 
win H-[z(A+1)] (4.100) 


to the left-hand side of (4.97), where H = [HT 0)" has been introduced. This reflects the 
special structure of the inverse dynamics problem. Applying feasible interpolatory quadra- 
ture formulas, leads to a computational form of (4.97) and (4.99). 


4.3. Numerical investigations 


This section aims to solve numerically the inverse dynamics problem of several mechani- 
cal systems, by applying the newly developed space-time integration methods. In partic- 
ular, the method of characteristics, termed MOC (Sec. 4.1), and the space-time finite ele- 
ment method, termed ST-FEM (Sec. 4.2), are applied. ST-FEM relies on bi-linear isopara- 
metric space-time elements. The evaluation of the integrals in (4.82) relies on standard 
Gaussian quadrature rules. This section partly reproduces [93]. 


Example 4.23 (Linear-elastic bar). Starting with the feedforward control of the linear- 
elastic bar, enables a comparison of the numerical results delivered by both MOC and ST- 
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FEM with the analytical solution (cf. Rem. 2.6 and Rem. 2.7). The prescribed trajectory of 
the bar at s = L is assumed as 


0 t < to 
1 R t — to T 1 

y(t) = 5 sin Ir zen - z) + 5 to <st<tr (4.101) 
1 t> tr 


for to = 1 and tf = 3. Moreover, the remaining data for this problem are EA = 1, pA = 1 
and L =1. 


MOC © 
analytic 


Figure 4.20.: Numerical solution for the actuating force f (t) computed using the method of characteristics 
(circles) compared with the analytical solution (solid line). 


ST-FEM ae STFEM © 
analytic analytic 


Figure 4.21.: Numerical solution for the actuating force f(t) computed using the space-time finite element 
method (circles) compared with the analytical solution (solid line). Left diagram: 5 x 15 = 75 elements, right 
diagram: 10 x 50 = 500 elements. 
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In Fig. 4.20 and Fig. 4.21 the numerical results for the actuating force f(t) are compared to 
the reference solution obtained by applying formula (2.22). It can be observed that the numer- 
ical results of the two alternative schemes under consideration closely match the reference 
solution. 


Figure 4.22.: Characteristic net along with a representative node Q located at (so, to) € Q. 


The results of MOC rely on a total of 356 unknowns. In Fig. 4.22, exemplarily a characteristc 
net is depicted. Note that due to the underlying hyperbolic structure, pre- and post-actuation 
phases have to be taken into account (cf: Rem. 4.9). Prescribing the motion on the boundary 
ats = L during the time interval [to, tf] leads to the requirement of pre- and post-actuation 
phases in the solution of the inverse dynamics problem. This can be observed from Fig. 4.22 
where the pre-actuation phase is related tot € [0, to] and the post-actuation phase is related 
tot € [t¢,T]. Since the slope of the characteristic lines depends on the velocity of wave 
propagation, the span of both, pre- and post-actuation phases is determined by the length 
of the bar along with the the velocity of wave propagation. The shaded trapezoidal area in 
Fig. 4.22 elucidates this correlation. 


Figure 4.23.: Functions u(s, t) = d;u(s, t) and n(s, t) = Bd,u(s, t) computed with the method of characteris- 
tics. Part of the characteristic net is also visualized in the s, t-plane. 
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Fig. 4.23 depicts the numerical solution for the two functions v(s,t) = Aır(s,t) and 
n(s,t) = Bdsr(s,t), respectively, computed with the method of characteristics. The corre- 
sponding characteristic net is also partially visualized in the s, t-plane. 


hı 


Figure 4.24.: Functions he (s, t) associated with the two Riemann Invariants. 


Figure 4.25.: Numerical solution for the displacement field u(s, t) and velocity field v(s, t) computed with the 
space-time finite element method. The corresponding finite element mesh (5 x 15 bilinear elements) can also 
be partially seen in the s, t-plane. 


Concerning ST-FEM, two alternative meshes comprised of rectangular bi-linear elements 
have been applied. The first one consists of 5 x 15 = 75 elements, leading to a total number 
of 195 unknowns. The second one consists of 10 x 50 = 500 elements, leading to a total 
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number of 1150 unknowns. Fig. 4.24 contains plots of the two functions ha(s, t) introduced 
in (4.47) associated with the Riemann Invariants (cf. Rem. 4.10). In particular, function 
hı(s,t) corresponds to the characteristic line Cı whose slope is given by ds/dt = k’(t) = c. 
Similarly, function hz(s,t) corresponds to the characteristic line C} whose slope is given by 
ds/dt = k' (t) = -c. It can be observed that both functions are indeed constant along the re- 
spective characteristic line. The numerical solution for the displacement field u(s, t) and the 
velocity field v(s, t) computed with the space-time finite element method are depicted in Fig. 
4.25. There, the corresponding finite element mesh is also partially visible in the s, t-plane. 


Concerning numerical accuracy, the method of characteristics yields exact results in the lin- 
ear setting. This is due to the fact that the difference scheme presented in Rem. 4.8 is capable 
to exactly integrate along straight characteristic lines. In contrast to this, the accuracy of 
ST-FEM is of order two. This can be concluded from Fig. 4.26, in which the relative error 


[pore = fre 
| 

is shown. Here, f°" is the actuating force calculated by applying the analytical solution 

(2.22), while f" refers to the numerical solution obtained with ST-FEM. To get the results in 


Fig. 4.26, the space-time finite element mesh has been uniformly refined. In this connection, 
he denotes the element length in space direction. 


ef (he) = 


107! 


1073 
107! 
he 


Figure 4.26.: Relative error in the actuating force f by using the space-time finite element method. 


Example 4.24 (Geometrically exact string - planar formulation). The next example deals 
with the planar motion of a geometrically exact string. In particular, the trajectory of the 
string at s = L is prescribed as straight line in the x, y-plane via 


y(t) = y(t) H (4.102) 
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where function y(t) is given by (4.101) with to = 2 and tf = 4. The remaining data is given 
by EA = 1, pA=1,L=1andg=9.81. 


The initial values for the inverse dynamics problem (2.11) have been obtained by solving 
the equilibrium problem of the string by means of the semi-discrete model described in Sec. 
3.1. For that purpose, 50 finite elements have been employed. The thus obtained vertical 
equilibrium configuration of the string is characterized by a total length (2.6) ofl(0) = 9.985 
and a suspension force of fọ = -g €y. 


Fig. 4.28 depicts the components of the actuating force required to realize the partially pre- 
scribed motion of the string. It can be observed that both the method of characteristics and 
the space-time finite element method yield closely related results. These results have been 
obtained by employing 50 x 150 = 7500 rectangular space-time elements leading to a to- 
tal number of 15600 unknowns. Merely 4 Newton iterations were necessary to achieve the 
present results by applying the simultaneous solution procedure outlined in Sec. 4.2. 


On the other hand, the results of the method of characteristics relies on a total number of 
6727 unknowns. Convergence was attained after 8 Newton iterations. It is worth noting that 
the data pertaining to the initial equilibrium configuration of the string has been used to 
initialize Newton’s method throughout the space-time domain, both for ST-FEM and MOC. 


t 


T 


0 L s 


Figure 4.27.: Characteristic net for the geometrically exact formulation of strings. 


In Fig. 4.27, exemplarily a charcateristc net is depicted. Note that similar to the linear setting, 
the solution of the inverse dynamics problem requires pre- and post-actuation phases. This 
is again indicated with the shaded area in Fig. 4.27. The parameters to, tr and T have the 
same meaning as in Ex. 4.23. 
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9.5 
0.25 


0.1 —9.7 
ak = 
9.9 
-0.1 
—0.25 -10.1 
0 2 : 4 6 0 2 ; 4 6 


Figure 4.28.: Actuating force components f£(t) and fy(t) at s = 0 accounting for a partly prescribed motion 
of the string at s = L. 


To get an impression of the overall motion of the string, snapshots in the x, y-plane and in 
space-time are shown in Fig. 4.29 and Fig. 4.30, respectively. In addition to that, the stretch 
distribution over the space-time domain is depicted in Fig. 4.31. 


10 + 


Figure 4.29.: Snapshots of the planar motion of the string in the x, y-plane. 
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Figure 4.30.: Snapshots of the motion of the string in the space-time domain. 


Figure 4.31.: Stretch distribution over the space-time domain. 
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Example 4.25 (Additional mass). In extension to Ex. 4.24, an additional point mass M = 1 
is attached to the right end of the string (cf. Rem. 2.2 as well as Rem. 3.2), then the term 
n(t) = M (g- ðv” (L, t)) has to be taken into account in weak form (4.80)2. Correspond- 


ingly, in the discrete formulation, the following entries of matrices (4.82), and (4.82); need 
be modified according to 


Aab — Aap + Ia mf Na: N, dt 
99, 


FE CFR +Mg | Na dt 
99 


for a,b € N}. The prescribed motion of the point mass is again governed by (4.102), with 
to = 1 and tr = 3. The remaining data is specified by EA = 10, pA = 1, L = 1 and g = 9.81. 


In a first step the equilibrium configuration of the string-mass system has been computed 
by applying 50 finite elements within the semi-discrete formulation (Sec. 3.1). Accordingly, 
the vertical equilibrium configuration providing the initial data for the subsequent inverse 
dynamics simulation is characterized by a total length (2.6) ofl(0) = 3.258 and a suspension 
force of fọ = —2g ey. 


The resulting components of the actuating force required to realize the partially specified 
motion are depicted in Fig. 4.32. 


Figure 4.32.: Actuating force components f,(t) and fy(t) at s = 0 accounting for a prescribed motion of a 
mass point attached to the string at s = L. 


It can be seen that the results of both methods under investigation match each other closely. 
The numerical results of the space-time finite element method are based on a mesh contain- 
ing 50 x 200 = 10000 elements leading to a total number of 20800 unknowns. Merely 4 
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Newton iterations are required to attain the converged solution. The corresponding results of 
the method of characteristics are based on a total number of 10515 unknowns. The required 
number of Newton iterations amounts to 7. Again the data pertaining to the initial equilib- 
rium configuration of the system has been used to initialize Newton’s method throughout 
space-time. 


The overall motion of the string is illustrated with snapshots of successive configurations of 
the system in Fig. 4.33. The total length (2.6) of the string versus time is depicted in Fig. 4.34. 
In addition to that, the stretch distribution over the space-time domain is depicted in Fig. 4.35. 


Figure 4.33.: Snapshots of the motion of the string with attached mass point. 
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Figure 4.34.: Total length of the string. 


110 


4.3. Numerical investigations 


Figure 4.35.: Stretch distribution of a string with attached mass point over the space-time domain. 


Example 4.26 (Geometrically exact string). This example aims to underline the impor- 
tance ofthe proposed method for the feedforward control of the three-dimensional motion of 
the string. In particular, a point mass, attached ats = L (cf. Rem. 2.2, Rem. 3.2 and Ex. 4.25) 
is supposed to move on a helix. To this end, the helix-shaped trajectory is prescribed by 


cos(27:p)-1 0 0 
y(t) =| sin@r-p) |Vtel y(t)=|0|Yt< to y(t)=| 0 |Wt>tp (4.103) 
zfp 0 Zf 


where 


1 1 1 t-t 1 
Ors cosl: y) +3 y= zain(x 2 -2)+ I= [totr] = [2,8] 


In Fig. 4.36 the prescribed coordinates are plotted versus time. The remaining data for this 
example is specified by EA = 10, pA = 1, L = 1, g = 9.81, M = 1 and zf = 5. 


In a first step the vertical equilibrium configuration of the system (Fig. 4.37) subject to gravity 
load is computed by applying the semi-discrete model of the string (cf. Sec. 3.1) consisting of 
15 finite elements. Accordingly, the vertical equilibrium configuration is characterized by a 
total length (2.6) of 1(0) = 3.258 and a suspension force of fọ = -2gez. The thus obtained 
equilibrium configuration provides the initial data for the inverse dynamics problem (2.11) 
to be solved subsequently. In this connection, the data is also used to initialize Newton’s 
method through the space-time domain. The space-time finite element approach (Sec. 4.2) is 
applied in conjunction with the simultaneous solution procedure outlined in Sec. 4.2. 
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Figure 4.37.: Equilibrium configuration of the string with attached point mass under gravity load. In addition 
to that, the helix-shaped trajectory to be tracked by the point mass is indicated as well. 


The computed three components of the actuating force required to realize the partially pre- 
scribed motion of the system are depicted in Fig. 4.38. Furthermore, Fig. 4.39 displays the 
total length of the string versus time. The resulting motion of the system is illustrated with 
Fig. 4.40 and Fig. 4.41. These results have been obtained by employing 15 x 149 = 2235 rect- 
angular space-time finite elements amounting to a total number of 7599 unknowns. Merely 
4 Newton iterations were required to reach the numerical solution. 


112 


4.3. Numerical investigations 


—17 


—18 


—21 


0 3 ó 9 0 3 + 6 9 


Figure 4.38.: Components of the actuating force f(t) obtained as solution of the inverse dynamics problem. 
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Figure 4.39.: Total length of the string. 
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-2 -1 0 
x 
Figure 4.40.: Computed trajectory of the actuated upper end of the string (s = 0) projected onto the horizontal 
x, y-plane. In addition to that, the projection of the prescribed trajectory of the lower end of the string (s = 1) 
is also shown. 


y -1 -2 


Figure 4.41.: Snapshots of the moving elastic string (the point mass attached at the lower end of the string is not 
shown). According to the inverse dynamics problem solved, the lower end of the string starts at rest (starting 
point y(to) = (0, 0, 0)) and traces the prescribed helix-shaped trajectory until the end point y(żf) = (0,0,5) 
has been reached. 


Example 4.27 (Recursive solution procedure). We next apply the recursive solution proce- 
dure described in Sec. 4.2. The data is kept the same as before. The number of finite elements 
in space direction is reduced from 15 to 10 in order to simplify the graphical illustration of 
the results. Accordingly, the space-time finite element mesh is comprised of 10x 149 rectan- 
gular elements leading to N = 10 time-space slabs used in the recursive solution procedure. 
As expected, the numerical results are practically indistinguishable from those obtained by 
applying the simultaneous solution procedure (cf. Sec. 4.2). In Fig. 4.42 the contact forces and 
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trajectories along the boundaries T, of the time-space slabs are shown. While the dimension 
of the iteration matrix pertaining to the simultaneous solution is equal to 5364, the dimension 
of the iteration matrix of the recursive solution procedure is equal to 1341 (cf. Rem. 4.20). On 
average 3.2 Newton iterations were required in each of the 10 steps of the recursive scheme. 
In comparison, 4 Newton iterations were required to reach the solution of the simultaneous 
solution scheme. 
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Figure 4.42.: Components of the contact-forces (left) and trajectories (right) along the boundaries of the time- 
space slabs obtained with the recursive solution procedure. 
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Remark 4.28 (Perturbed data of the prescribed trajectory). According to Hadamard, a 
problem is well-posed if there exists a solution to the problem which is unique and stable 
(see [45] or [58, Def. 1.13]). Since existence and uniqueness of the solution to the present 
space-time boundary value problem are assured by the differential flatness of the system, 
stability still has to be ensured. The solution is called stable if it depends continuously on the 
given data - small perturbations of the prescribed servo-constraints on Ny must not lead to 
unbounded results for the actuating force f(t) on 9Qy. This can be verified numerically by 
adding random noise to the given data of the trajectory such that 


’-»l,=® 


where yê denotes the perturbed trajectory. The perturbed components of y° are shown in Fig. 
4.43. For comparison, the original components of y given by (4.103) are depicted in Fig. 4.36. 
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Figure 4.43.: Components of the perturbed data for the prescribed trajectory yê for ô : t œ [0, 0.05]. 
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Fig. 4.44 shows the results for the actuating forces resulting from the solution of the inverse 
problem with the perturbed data. Comparison with Fig. 4.38 indicates that the perturbed 
problem still yields results close to the original ones. 


a=0.02 --- a=0.02 --- 
a=0.05 — a=0.05 — 


a=0.02 --- 
a=0.05 — 


Figure 4.44.: Components of the actuating force obtained as solution of the inverse problem with noisy data 
for the prescribed trajectory yê at 90, with ô : t > [0,a] anda € {0.02, 0.05}. 


Example 4.29 (Geometrically exact beam formulation - circular trajectory). A beam, 
with mass density p = 1 and axial-, bending- and shear stiffness EA = 1, EI = 1 and GA = 1 
respectively, is investigated. The length of the beam in a stress-free reference configuration 


is assumed to be L = 1. The aim is, to find the actuation f = Ee fy m|" applied ats = 0, 
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4. Simultaneous space-time integration 


such that the end-effector of the beam ats = L follows a prescribed trajectory. For this, a 
circular motion starting at tọ = 2 and ending at tf = tọ +T = 4 is choosen as given by 


1 1-cos( - y(t)) 0 
y= |0| Vt < to y=| sin(g-y(t)) | Vt € [to to+T] y= |1] Vt > tp. (4.104) 
0 o- y(t) p 


Herein p € [0,22] has been introduced, specifying the angle of rotation. Furthermore, the 
function y : T > R is obtained as 


0.6 
0.2 
= = 
-0.2 
-0.6 
-0.5 
0 2, 4 6 0 2, 4 6 


Figure 4.45.: Components of the force (left) and torque (right) acting at s = 0 computed with the proposed 
space-time finite element method such that the beam at s = L follows a prescribed circle from Po(1,0) to 
Pr(0, 1). 


Figure 4.46.: Snapshots of the motion of a geometrically exact beam actuated at s = 0 such that its end at s = 1 
aligns with a prescribed circular motion. 
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4.3. Numerical investigations 


In Figure 4.45 the components of the actuating force fi (left) and the actuating torque m 
(right) is depicted. Note that also the delay time can be observed herein. This is due to the 
hyperbolic structure of the underlying system mentioned earlier. In Figure 4.46 snapshots of 
the planar motion of the beam satisfying the prescribed trajectory ats = L are shown. These 
results have been obtained by employing 10 x 50 rectangular space-time elements leading to 
a total number of 3450 unknowns. 


Example 4.30 (Geometrically exact beam formulation - straight trajectory). Considering 
a beam with material and geometrical properties taken from Ex. 4.29, is forced ats = 0 such 
that the beam at s=1 follows the prescribed motion 


1 1+ y(t) 2 
y = |0| Vt < to y=} y(t) | Vt € [to to+T] y= |1|Vt > tp. (4.105) 
0 0 0 


where the function y : T +> R is assumed as 


y(t) =1- ; (cos (= (sin (x- — _ =) + 1)) + 1) 3 


A simultaneous space-time integration of the partial differential equation at hand leads to 
the actuating force and torque at s = 0, searched for. 


0.2 


0.5 


Silt) 
m(t) 


—0.5 —0.2 


t t 


Figure 4.47.: Components of the force (left) and moment (right) acting at s = 0 computed with the proposed 
space-time finite element method such that the beam at s = L follows a prescribed straight line from Po (1, 0) 
to Pr (2,1). 


In Fig. 4.47, the force components fy and fy as well as the torque m are depicted. These results 
have been obtained by employing 10 x 50 rectangular space-time elements leading to a total 
number of 3450 unknowns. 
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4. Simultaneous space-time integration 


0 1 2 


Figure 4.48.: Snapshots of the motion of a geometrically exact beam actuated ats = 0 such that its end at s = 1 
aligns with a prescribed straight motion. 


Example 4.31 (Nonlinear continuum - plane stress). This example aims to demonstrate 
the relevance of the newly proposed simultaneous space-time integration of the inverse dy- 
namics problem. For this, a quadratic membrane having edge length l = 0.5, mass density 
p = 1, Young’s Modulus E = 100 and Poissons ratio v = 0.3 is considered. The actuation 
f :3Qf XT + R” is searched for, such that y : OQ, x T > R” can be prescribed as 


1 1+ y(t) 2 
y = |s2| Vt < to y=} y(t) | Vt e [t,t +T] y=|s2+1[Vt>tr. (4.106) 
0 0 0 


In Eq. (5.13), the function y : T > R is given by 


y(t) =1- > (cos (= (sin (x- — - =) + 1)) + 1) 


Note that the fully boundary-actuated problem, that has been considered in Ex. 4.23, Ex. 
4.24, Ex. 4.25, Ex. 4.26, Ex. 4.29 and Ex. 4.30, is now extended to partly boundary-actuated 
problems, e.g. 057 U 05, # ƏS is explicitly permitted. 
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4.3. Numerical investigations 


Figure 4.49.: Inverse computation of the actuating traction on dQ ¢ accounting for a partly prescribed motion 
on dQy. 
n 


In Fig. 4.49 the components of the actuating force fi : Qf X T > R, computed using the 
space-time finite element method, are depicted. In Fig. 4.50 snapshots of the prescribed ma- 
neuver are illustrated. These results have been obtained by employing 10x10x80 hexahedral 
space-time elements leading to a total number of 40480 unknowns. 


t=0.00s 


Figure 4.50.: Snapshots of the membrane actuated on 90; such that the boundary 9Q,, translates by 


u= [1 1 0] ‘3 by simultaneously keeping its shape straight. 
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5. Flexible multibody systems 


After the inverse dynamics problem for flexible mechanical systems, including slender 
structures such as strings and beams but also general three dimensional continua under- 
going large deformations, has been analyzed, simultaneous space-time integration tech- 
niques could have been established in order to solve the control problem at hand. In the 
present Chapter, we aim to apply the newly established theory to the feed-forward control 
of flexible multibody systems. Therefore, the newly proposed simultaneous space-time 
integration is first adopted in Sec. 5.1 to the cooperative control problem of rigid bodies 
actuated through several flexible strings. Later on in Sec. 5.2, a rotational arm model 
consisting of two rigid and one flexible arm is investigated. Accompanying numerical 
examples are underlining the relevance of the proposed methods in context of flexible 
multibody systems. 


5.1. Cooperative control problem 


This section is concerned with finding the feed-forward control of rigid bodies actuated 
through multiple elastic strings. More precisely, four elastics strings are attached to a 
rigid body. The task is to find the actuating forces that are applied to the free end of each 
rope such that the rigid body follows a prescribed motion. To avoid singularities, we 
presume that the contact points do not lie on a plane, a straight line or are concentrated 
in one point. Instead we assume the contact points to be vertices of a tetrahedron. See 
Fig. 5.1 for an illustration of the control task. 


The structure of the equations of motion of the underlying multibody system, comprised 
of geometrically exact ropes and rigid bodies, is outlined first. Numerical examples will 
be given underlining the applicability of the proposed solution strategy. 
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4 


Figure 5.1.: Elastic string (left) as sub-system of a multibody system (right): Cooperative transportation of a 
rigid body through four elastic strings actuated by external forces fp, k € {1,..., 4}. In the inverse dynamics 
problem, the motion of the rigid body is prescribed and the goal is to determine the motion of the whole 
multibody system along with the actuating forces f. The prescribed motion of the rigid body gives rise to 
the boundary conditions (trajectory and contact force nX) for each string at its point of attachment to the rigid 
body. 


Rigid body. When considering general rigid bodies, in principle the same strategy as 
for the attached mass point can be applied (cf. Rem. 2.2 and Ex. 4.24 as well as Rem. 2.3). 
The actuating forces needed to achieve the desired motion of the rigid body can then be 
calculated directly by considering the governing equations of motion. For that purpose, 
the rigid body is modelled as a Cosserat continuum subjected to geometric constraints. 
Consequently, the motion of the rigid body is fully described by the directors d; : T > R? 
and the position of the centre of gravity x : T +» R?. Furthermore, the rigid body 
is assumed to have total mass M = te po dV with a density function pọ : Bo > R. 
Accordingly, the equations of motion can be written as 
Mé?x — fxs -G=0, 
Eijörd; = ‘Picts + Ajjdj =0, (5.1) 


gc(di) = di -dj — 6); =0. 


Herein G = Mge; is the gravitational force and E;; are the components of the referential 
Euler tensor which is closely related to the classical inertia tensor of rigid body dynamics 
(see [21] for more details). Furthermore, f,,., is the resultant external force and 


f= xine (5.2) 


are the external director forces (Fig. 5.2). This holds for a discrete actuation of the rigid 
body. Note that by introducing the matrix 


I 
X*I 
xt)’ 
XkI 


Hjk = (5.3) 
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5.1. Cooperative control problem 


containing the information of the geometric position of the contact points, the follow- 
SR: f j r AT 
ing linear relation between the external (director) forces f(t) = [ Fic arl and the 


contact forces nk(t) is obtained as 
F(t) = An). (5.4) 


The geometric constraints (5.1)3 are enforced by the Lagrange multipliers A;;. Addition- 
ally to the standard holonomic constraints (5.1)3, control constraints as given by 


9;=4- y(t) =0 (5.5) 


are imposed to force the rigid body at hand to follow a prescribed motion y(t) = [yx Ya ] R 
In Eq. (5.5), the function g : T > R!? has been defined as q = [x di] .T The servo- 
constraints (5.5) of course must not violate the holonomic constraints (5.1)3. The differ- 
ential part of system (5.1) along with the control constraint (5.5) gives rise to a purely 
algebraic system of equations given by 


Flt) = Dažy(t) + Fy(t) -6. (5.6) 
Herein, the coefficients D, F and G has been defined as 
M 0 0 0 = G 
alerts on 


Eq. (5.6) determines the external (director-) force f(t) that is conjugate to the prescribed 
trajectory y. Once f(t) has been computed, the k € N contact forces n*(t) for the k 
strings at sk = 1 can be easily computed by knowing the position of the contact point of 
the string at the rigid body. By defining the components Qx; of the inverse of the positive 
definite matrix H € R441 introduced in (5.3), a linear relation of the contact forces to 
the (director-) forces is obtained as 


nk (t) = Qu f(t). (5.8) 


Once the forces n*(t) have been calculated, each string can be solved separately by in- 
serting the forces into the corresponding Neumann boundary condition (1.6) of the initial 
boundary value problem established in Sec. 1.1. This cascade-like solution strategy for the 
cooperative control problem enables a parallelisation of the computation of the inverse 
dynamics of the attached strings. The applicability of this approach is demonstrated in 
Ex. 5.3. 
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5. Flexible multibody systems 


Figure 5.2.: Force nx acting on a rigid body at point P = x + X;dj. 


Remark 5.1 (Lagrange multipliers). The actuating forces depend on the Lagrange multi- 
pliers Aij, i.e. an unique solution for the actuation of the rigid body requires 6 independent 
components Àk : T +> R fork € {1,--- ,6} C N of Aj; to be partly specified. 


Ak(t) = ya, (t) 


Essentially this amounts to partly specifying the stresses within the rigid body (cf. [55]). 


Multibody system. For the purpose of verification, the flexible multibody system is 
introduced for the direct dynamics problem. Therefore, consulting (3.4) along with (3.5) 
and (3.7) for the string and (5.1) for the rigid body, the motion of the overall system is 
governed by a differential-algebraic system of equations given by 


0 = R(D7q. 9, t) - (349; (9))" u 
0= g; = V p+ _ (x +X;d;) 
0 = ge = di -dj — ôij 


where 


2 z n |MD2q+fi®-G 
RDGGAN=| ppg+rg-G 


Therein the solution q is defined by q = Ci 7] " Note that the Lagrange multipliers 


u: T > R? in conjunction with the constraint Jacobian 9g9,(g) can be identified as the 
contact force acting among the string and the rigid body 


T 
pðag() = [0 >> p -u -4Xı -4X2 -uXs] (5.9) 


See Ex. 5.2 for an application of this model to a simulation of a free swinging rigid body 
attached to a single elastic string. 
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5.1. Cooperative control problem 


Example 5.2 (Direct dynamics problem - pendulum). To verify the numerical implemen- 
tation of the procedure outlined in this section, an example for the forward dynamics of a 
rigid body suspended by a flexible string is considered. Therefore, a cube with edge length 
a = 2 and total mass M = 0.1 is attached to an elastic string characterized by its Young’s 
Modulus E = 1, density po = 1, cross-sectional area A = 1 and length L = 1. The velocity 


v(t) of the rigid cube at time t = 0.0 s is assumed to be vp = [2 -5 o|’. In Fig. 5.3 
snapshots of the swinging rigid body are depicted. 


Figure 5.3.: Free oscillation of a rigid cube with edge lenth a = 2 and total mass M = 0.1 connected to a flexible 
rope with p = 1 and E = 2. 


Example 5.3 (Inverse dynamics problem - rotation and translation). To verify the pre- 
sented approach to the inverse dynamics problem under consideration, the following scenario 
is given. A rigid cube with edge length a = 2 and mass M = 1 is supposed to accomplish 
a rest-to-rest maneuver. For this a translation of the cube from point Py = (0,0,0) to point 
Pr = (2,2,2) along a straight line together with a simultaneous rotation of x around the 
z-axis is planned. The maneuver is intended to start at tọ = 1.0 sand end at tr = 9.0 s. 
The motion of the rigid body described in terms of a Cosserat point subjected to geometric 
constraints can then be prescribed together with 


0=-- [cos (= (sind +1)) - 1] for Zn -;| (5.10) 
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5. Flexible multibody systems 


through the trajectory as given by 


-cosp+1 
-cosp+1 
-cosp+1 
cos Q -1 
sin o 0 


N N 


Vt<to Y= ; VteeTm=ltiotrl, Y, = Vi>tr 


Yro = 


Hooo-.o0o00r 000 
oO 


(5.11) 


Thereby, four elastic ropes with mass density p = 1 and Young’s Modulus E = 1, are used to 
actuate the rigid body for which the motion is prescribed. Regarding the k predefined contact 
points fori € {1, 2, 3} 


P =x+X*d; 


and following Sec. 5.1 the contact forces needed to achieve the desired motion can then be 
computed using (5.8) with 


I I I I 

I -I I -I 
Fr I -I -I I 
-I I I I 


where I € R%@ is the identity matrix. Note that, following Rem. 5.1, the Lagrange multipliers 
cannot be left undefined. For the given maneuver, a uniaxial tension within the rigid body 
is choosen: 


A=[0 0 Mg 0 0 of (5.12) 


In Fig. 5.4 components of the k contact forces are shown. The forces acting on the upper 
end of the k ropes at sk = 0 can be calculated such that the rigid body at hand follows the 
prescribed trajectory. The numerical solution is shown in Fig. 5.5. 


To verify the outlined method, the forces f,(t) can be inserted into the flexible multibody 
system presented in Sec. 5.1. In Fig. 5.6, snapshots of a forward simulation of the flexible 
multibody system at hand actuated with the forces f(t) acting at the upper ends of the k 
ropes computed numerically using the approach presented above are shown. 
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5.1. Cooperative control problem 


Figure 5.4.: Components of the contact forces nx (t) such that the rigid body follows the prescribed translation 
of the rigid body from point Po = (0, 0,0) to point Pr = (2, 2,2) accompanied by a simultaneous rotation 7r. 
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Figure 5.5.: Components of the actuating forces f (t) for k € {1, 2, 3,4} acting on the upper end of the ropes 
at sk = 0 accounting for the prescribed motion of the rigid body. 
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5.2. Rotational arm model 


t=10.0s 
t=7.0s 
Hs t=6.0s 
t=0.0s 


Figure 5.6.: Snapshots of a rigid cube with edge length a = 2 and total mass M = 1 actuated through four 
flexible ropes with p = 1 and E = 2 following a straight line from Po = (0, 0,0) to Py = (2, 2,2) while rotating 
simultaneously by the prescribed angle z. 


5.2. Rotational arm model 


A rotational arm model consisting of a chain of three arms is investigated subsequently. 
Particularly, a flexible arm 63 is attached to a rigid arm B; in point C, which in turn is 
connected to another rigid arm B, in point B. Body $; in turn is attached to an inertial 
frame in point A. See Fig. 5.7 for an illustration thereof. The bodies are assumed to have 
density p; : B; > R and total mass M; = Jo, pi dV for i € {1,2,3}. The task is to find 
the torques m4, mg and mc that need to be applied to the joints A, B and C, respectively, 
such that the end-effector P follows a prescribed motion y(t) given by 


1 1 — cos(© - y(t)) 0 
y = |0| Vt < to y=] sin(©-y(t)) | Vt € [to,to+T] y= |1| Ve > tr . (5.13) 
0 ©- y(t) © 


In Eq. (5.13) the angle of the prescribed overall rotation © € [0,27] and the function 


y(t) =1- - (cos (= (sin (x- = - =) + 1)) + 1) 


has been introduced. 
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Figure 5.7.: Rotational arm model. Geometrically exact beam in context of flexible multibody system (2d). 
The coordinates (xp, yp) as well as angle Op of the end-effector P denote the system outputs. The torques 
(ma, mg, mc) in joints (A, B, C) are the conjugate system inputs. 


In principal, the same cascade-like procedure, that could have been developed in Sec. 5.1 
in context of the cooperative control problem, can be applied to the inverse dynamics of 
the flexible rotational arm model at hand. The strategy reads as follows. 


The flexible arm is modelled as a Cosserat rod aligning with the framework introduced in 
Sec. 1.1 and Sec. 2.2, analyzed in Ex. 4.5 as well as Ex. 4.4 and solved numerically in Ex. 
4.30 and Ex. 4.29. The proposed numerical schemes based on a simultaneous space-time 
integration, do not only provide a solution for the actuating force searched for. They also 
give rise to the overall motion, including the trajectory of the actuating joint C as given 
by 

xc(t) 
yc(t) 
In case of describing the kinematics of the rigid part of the system in Cartesian coordi- 
nates, i.e. in a redundant configuration space, the trajectory yç can directly be inserted 
into the differential algebraic system of equations. Similarly to Sec. 5.1, the external 
applied (director) forces need then be determined. 


Yc) = ; (5.14) 


Alternatively, minimal coordinates can be used to describe the kinematics of the rigid 
multibody system. Then we need to find coordinates ©; and Og, aligning with the com- 
puted trajectory of the sub-systems end-effector (joint C), i.e. 


xc(t)} _ | cos O(t) + cos ©2(t) 
yc(t)| — | sin @,(t) + sin ©2(t) 


Obviously, the inverse kinematics (5.15) are not uniquely solvable! for the minimal coor- 
dinates ©, (t) and ©2(t). We therefore solve (5.15) in terms of a least-square ansatz: 


y(t) = = F(0,,9) . (5.15) 


|F @, Q2) - Yell — Min. 


1 By considering a feasible configuration ©; = ©, and ©, = ©; such that 7(t) = y(t) but x(t) # xı (t) and 
X2(t) # x2 (t) proves the singularity of (5.15) 
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Once the motion ofthe remaining rigid part ofthe model has been calculated, the actuat- 
ing torques, can be determined uniquely by evaluating the governing balance equations 
of linear and angular momenta of the rigid part ofthe model. 


Since the contact force nc acting at joint C is determined by the inverse computation of 
the flexible beam 83, the contact force in joint B can be obtained, by using the balance of 
linear momentum of the rigid body 62, as 


ng = D;P; = G, —-nc. (5.16) 


In Eq. (5.16), the linear momentum of body Bz has been introduced as D;P2 = M2a2. 
Once the contact force ng is determined, the balance of angular momentum of the rigid 
body B, provides the contact torque mg acting at joint B 


Mc — mg + (r2 X G2) + (rg x ng) + (rc X Nc) = D;La, . (5.17) 


Then the torque my acting at joint A can be calculated by establishing the balance of 
angular momentum of body 5; 


ma + mp + (rı X G1) + (rg X Np) =D;La, . (5.18) 


In Eq. (5.17) and Eq. (5.18), the gravitational forces G; and Gz, acting at the center of 
gravity of body Bı and By, respectively, have been introduced. The torques m4 and mg 
applied to joints A and B, respectively, can be identified as differentially flat inputs cor- 
responding to the differentially fat output, by means of the trajectory of the end-effector 
of the fully actuated rigid multibody system. In Fig. 5.8 the inverse computed actuating 
torques m4, mg and mc are depicted. 


Figure 5.8.: Actuating torques m; for i € {A,B,C}, applied to joints A, B and C such that the end-effector 
captures the prescribed motion. 
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6. Conclusion 


6.1. Summary 


The present work has focused on the inverse dynamics of underactuated flexible mechani- 
cal systems undergoing large deformations. Therefore, it could be demonstrated that this 
specific problem is governed by quasi-linear second-order partial differential equations 
subjected to time-varying Dirichlet boundary conditions that are enforced by spatially 
disjunct Neumann boundary conditions. Various mechanical systems including slender 
structures such as strings and beams as well as general three-dimensional continua have 
been presented and proved to align with this specific initial boundary value problem. 


Then, a brief survey on common solution strategies relying on a sequential space-time 
discretization has been given. Therefore, it could be elaborated that such methods are 
leading either to (unstable) internal dynamics or high differentiation index’ along with 
high demands on the smoothness of the prescribed constraint function. Both hinder a 
numerical stable integration of the semi-discrete equations. We were able to demonstrate, 
that these issues are caused by the non-standard construction ofthe constraint realization 
in conjunction with the applied sequential space-time discretization. 


At this point, an in-depth analysis of the initial boundary value problem at hand was un- 
dertaken. In that process, characteristic manifolds along which information flows could 
be identified and the hyperbolic structure of the governing partial differential equations 
has been revealed. Enlightning these wave phenomena has provided precious insights 
into the governing equations and has paved the way to develop novel numerical meth- 
ods that are capable of solving the inverse dynamics of (infinite-) dimensional systems. 


Therefore, a method that is based on a geometrical interpretation of the underlying par- 
tial differential equations could be introduced. It has been demonstrated concisely that 
once characteristic manifolds have been found, the partial differential equations at hand 
can be transformed equivalently into a system of ordinary differential equations along 
these manifolds. The succesful implementation of this method has motivated to develop 
a novel Galerkin-based space-time finite element method. This space-time finite element 
approach relies on continuous test functions. This is in sharp contrast to previously de- 
veloped space-time finite element methods that rely on discontinuous test functions thus 
making possible to divide the space-time domain into a sequence of space-time slabs even- 
tually recovering the time-marching format commonly applied in structural dynamics. 
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Then, these novel numerical methods could be adopted to the feed-forward control prob- 
lem of flexible multibody systems. Therefore, a rigid body that is cooperatively controlled 
through several elastic strings and a rotational arm consisting of rigid and flexible arms 
has been considered. 


Selected representative numerical examples have been presented demonstrating the ca- 
pabilities of the newly devised space-time integration methods. 


6.2. Outlook 


This work has answered some important questions but also raised new ones. Some of 
them are briefly summarized below. 


The simultaneous space-time integration methods proved to be the key to succesively 
solve the inverse dynamics problem of infinite-dimensional systems. Thereby it turned 
out, that the method of characteristics is very well adapted to the governing equations. 
On the other hand the space-time finite element method has proved to be versatile appli- 
cable due to its inherent simplicity. However, its potential has not been fully exploited yet. 
For instance, the possibility of using higher order elements as well as using unstructured 
meshes along with an error-based adaptive mesh refinement have not been implemented 
so far. Moreover a combination of the method of characteristics and the space-time finite 
element method might be conceivable by orienting the finite elements along characteris- 
tic lines. Such a characteristic Galerkin method may combine the accuracy of the method 
of characteristics and the simplicity of the space-time finite element method. 


Besides the spatially continuous systems, spatially finite-dimensional systems constitute 
another important class of inverse dynamics problems. Such systems are governed inher- 
ently by ordinary differential equations and hence elude from the framework proposed 
for spatially infinite-dimensional systems. Since the inverse dynamics of spatially dis- 
crete sytems such as cranes or manipulators with passive joints are highly relevant, two 
global Galerkin-based methods could have been developed in Rem. 4.21 and Rem. 4.22. 
One is derived directly from the continuous space-time finite element method presented 
in Sec. 4.2 and the other one is based on discontinuous approximations of the test- and 
trial functions along with additional flux terms propagating backwards in time. These 
fluxes are accounting for the causality of time-dependent physical systems. Further anal- 
ysis seams to be promising. 


Throughout this work, the motion of mechanical systems has partly been specified. But 
for some cases it might be beneficial to find trajectories connecting two prescibed points 
that are optimal in terms of causing a minimal effort of the actuation. Questions on the 
existence and uniqueness of a trajectory quickly arise, hindering a stable numerical solu- 
tion. Therefore, adopting the newly gained insights, regarding the property of physical 
time-dependent delay systems and its connection to the property of differential flatness 
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of spatially discrete systems, to the theory of optimal control, might be convenient to find 
feasible trajectories. 
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